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Abstract 

The analysis of the electronic properties of strained or lattice deformed graphene combines ideas from 
classical condensed matter physics, soft matter, and geometrical aspects of quantum field theory (QFT) in 
curved spaces. Recent theoretical and experimental work shows the influence of strains in many properties 
of graphene not considered before, such as electronic transport, spin-orbit coupling, the formation of Moire 
patterns and optics. There is also significant evidence of anharmonic effects, which can modify the structural 
properties of graphene. These phenomena are not restricted to graphene, and they are being intensively 
studied in other two dimensional materials, such as the transition metal dichalcogenides. We review here 
recent developments related to the role of strains in the structural and electronic properties of graphene and 
other two dimensional compounds. 
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1. Introduction 

Graphene has a number of special properties [1]. The unique combination of a two dimensional (2D) 
structure and massless electronic carriers leads to a variety of phenomena never considered before. These 
features are amplified by the high stiffness of the graphene lattice. The analysis of the electronic properties 
of strained or lattice deformed graphene combines ideas from classical condensed matter physics, soft matter 
[2, 3, 4], and geometrical aspects of quantum field theory (QFT) in curved spaces [5]. 

A vast literature already exists on topics such as ripples, effective gauge fields and ’’strain engineering” 
in graphene [6] and there are already a number of excellent books and reviews describing the advances in 
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the understanding of strain related effects [7, 8, 9]. The understanding of these phenomena has motivated 
the study of related effects in other 2D materials. 

Recent theoretical and experimental work shows the influence of strains on many properties of graphene 
not considered before, such as electronic transport, spin-orbit coupling, the formation of Moire patterns and 
optics. There is also significant evidence of anharnronic effects, which can modify the structural properties 
of graphene. These phenomena are not restricted to graphene, and they are being intensively studied in 
other 2D materials, such as the metallic dichalcogenides. 

We review here recent developments related to the role of strains in the structural and electronic prop¬ 
erties of graphene and other 2D compounds. We emphasize how strains influence a number of features not 
considered before. 

The review is divided into five self-contained sections which cover the main recent developments. A 
summary of basic properties, formalism and notation used throughout the review are found in Section 2. 
This section introduces the formalism by which elastic deformations of the lattice encoded in the strain 
tensor of the elasticity theory couple to the electronic degrees of freedom described by the Dirac equation in 
the continuum limit. A symmetry approach is used to derive all possible electron-phonon couplings and their 
physical content is described. The ’’elastic” gauge fields and the deformation potential to be used through 
the review are defined there. This section also includes a geometrical description based on quantum field 
theory in curved spaces that allows to predict some new effects as the Fermi velocity as a space-dependent 
tensor. This formalism is also used to describe some new developments on topological properties of the 
materials. 

Section 3 is devoted to the new developments concerning the membrane properties of the lattices. Al¬ 
though these aspects attracted attention from the birth of graphene [9], recent experiments [10, 11, 12] 
have shown that anharmonic effects in graphene and other 2D crystalline membranes have been seriously 
underestimated. This section describes the present situation concerning the effects of these anharmonicities 
on the elastic and mechanical properties of the membranes. A summary of the classical theory of flat mem¬ 
branes (Section 3.1.1) to fix the notations and definitions, is followed by a description of the modifications 
to the electron-phonon interaction due to defects (Section 3.1.2) and quantum corrections (Section 3.1.3). 
The modifications to the mechanical and thermodynamic properties are described in Sections 3.2, 3.3. The 
section concludes with an overview of the recent advances on the influence of topological defects on the 
membrane properties (Section 3.4). 

Section 4 deals with the new developments concerning the effect of strain on the electronic properties 
of the samples. Particular attention is put on the analysis of the correlation between carrier mobilities 
and amplitudes of strain distributions demonstrated recently (Sec. 4.1). The combined influence of strain 
and interactions is analysed in Sec. 4.5 which addresses in particular the access to topological phases 
through the interplay between strain and electron-electron interactions. Interaction induced magnetic or 
superconducting phenomena in strained graphene are also described. The optical properties of graphene and 
other atomically thin materials promise a wealth of applications. In particular, graphene plasmonics in the 
THz and mid-infrared bands [13], as well as graphene based broadband devices working in the near-infrared 
and visible part of the spectrum [14], have recently attracted attention. Further, the unusual magneto¬ 
optical properties of graphene, which are determined by the nonlinear Landau level spectrum of massless 
Dirac Fermions [15, 16], offer a great potential for the development of non-reciprocal optical devices. Optical 
properties are described in Sec. 4.6. 

Isolated atomic planes can now be reassembled into heterostructures made in a precisely chosen sequence 
[17]. Section 5 deals with the electronic properties of these new compounds that include superlattices of 
graphene deposited on other 2D crystals like hBN, SiC or graphene itself. The induced strains on the pristine 
sheet, resulting from a competition between the adhesion potential and the elasticity of the layer, are the 
factor ultimately determining the electronic and mechanical properties of the sheet. A general description 
(5.1) is followed by an analysis of the spontaneous deformations arising in the samples and their experimental 
consequences (5.2). 

Section 6 is devoted to the new 2D materials whose synthesis followed naturally that of graphene. It 
includes a bibliographic summary of the on-going research in isolated mono- and few-layers of hexagonal 
boron nitride (hBN), molybdenum disulphide, other dichalcogenides, black phosphorus and layered oxides. 
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There are important aspects of the vast subject of physics of strain that have been left out of the review. 
Some of them have been appropriately described in previous literature [8]; others, including mechanical and 
numerical approaches to analyse the stability of the lattice [18], have been omitted for lack of expertise. 
This review is the collective work of a theory group in Madrid that has been active in the field from the 
beginning and continues working on it. 

2. Geometric aspects: Group symmetry approach 

The aim of this chapter is to discuss the minimal models for the low energy modes of graphene. Both 
symmetry group and differential geometry arguments are employed in order to construct effective actions 
for the dynamics of graphene electrons and phonon modes, and the coupling between them. The standard 
k ■ p derivation is compared with a quantum field theory approach based on the study of the Dirac equation 
in curved spaces. Finally, we review some topological aspects of graphene physics, with special emphasis on 
topological defects, closely related to the structural properties of the graphene lattice. 

2.1. Electrons and phonons in graphene. Symmetry approach 

The idea of constructing effective actions for physical systems based only on symmetry considerations 
has been extensively used both in quantum field theory (QFT) and in condensed matter physics and it lies 
at the heart of the Landau Fermi liquid theory of metals. In solid state systems, it is the underlying crystal 
group what imposes the symmetry constraints on the action describing the dynamics of both electrons and 
phonons. 

Graphene consists on a single layer of sp 2 hybridized carbon atoms forming a honeycomb lattice, which 
is a triangular Bravais lattice with two atoms per unit cell. Three of the four electrons in the outer-shell 
participate in the strong a bond which keeps them covalently attached forming a planar structure with a 
distance between atoms of acc = 1-42 A. These a electrons are responsible for the structural properties of 
graphene, in particular its stiffness. The remaining electrons occupying the p z orbital perpendicular to the 
graphene plane are free to hop between neighbouring sites, leading to the ir bands. In pristine graphene, 
the Fermi level lies at the two inequivalent corners of the hexagonal Brillouin zone, K± (see Fig. 1). 

The three dimensional (3D) point group of the graphene crystal is D 6 h = D 6 ® i = Cq v ® ov We follow 
the notation of Ref. [19], which can be simplified if we only refer to the operations of the planar group, Cq v , 
that contains 12 elements: the identity, five rotations along the axis perpendicular to the graphene plane 
(z axis) and six reflections in planes perpendicular to it. Additionally, we have to take into account the 
reflection operation along the 2 axis, erg, when dealing with out-of-plane distortions. 

Instead of dealing with degenerate states at two inequivalent points of the Brillouin zone one can enlarge 
the unit cell in order to contain six atoms. Therefore, the folded Brillouin zone is three times smaller and 
the K± points are mapped onto the T point, see Fig. 1. From the point of view of the lattice symmetries, 
this means that the two elementary translations (tg 1 , tg 2 ) are factorized out of the translation group and 
added to the point group, which becomes C' ( ' v = C 6v +tg 1 x C 6v + tg 2 x C 6v [20]. This allows us to treat 
all electronic excitations at the corners of the Brillouin zone on a same footing, including also possible 
inter-valley couplings. The Bloch wave function is given by a 6-component vector which represents the 
amplitude of the p z orbitals at the 6 atoms of the unit cell. This vector can be reduced as A\ + B 2 + G'. The 
1-dimensional irreducible representations A\ and B 2 correspond to the bonding and anti-bonding states at 
the original T point, whereas ip ~ G' corresponds to the Bloch states at the original Brillouin zone corners. 
Then, in order to construct the electronic Hamiltonian for quasiparticles around the K± points, we must 
consider the 16 hermitian operators acting in the 4-dimensional space defined by the Bloch wave functions 
ip. These operators may be classified according to the transformation rules under the symmetry operations 
of Cq v , taking into account the algebraical reduction 

G' xG' ~A 1 + A 2 + B 1 + B 2 + E 1 + E 2 + E[ + E' 2 + G’. 

In order to make the discussion more clear we introduce the basis ip = (ipA+,ipB+,4’A-,fpB-), where each 
entry represents the projection of the Bloch wave function around each valley I\± on sublattice A/B. Then, 
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Figure 1: Graphene Brillouin zone corresponding to a real space unit cell with with 2 (dashed line) and 6 (solid line) atoms in 
the unit cell. 
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Table 1: Classification of the electronic operators according to the irreducible representations of Cq v . The signs zb denote if 
the operator is even or odd under time reversal. 


we introduce two inter-commutating Pauli algebras cr, and r, ; associated to sublattice and valley degrees of 
freedom respectively. The 16 possible electronic operators are generated by considering the direct products of 
the elements of these algebras (and the identity). Their symmetry properties are summarized in Tab. 1. We 
must take into account also the time reversal operation, which is implemented by the antiunitary operator 
T = t x IC, where K, represents complex conjugation. 

The effective low energy Hamiltonian is constructed as an expansion in powers of the crystalline momen¬ 
tum k = ( k x , k y ) ~ Ei around K± points. Up to first order in k we have 

H 0 = v F ( r z a x k x + a y k y ), (1) 

where v F is the Fermi velocity. This Dirac Hamiltonian describes the approximately conical dispersion of 7r 
bands around K± points. 

Phonon modes of graphene within the tripled unit cell can be also classified according to the irreducible 
representations of Cq v as shown in Tabs. 2 and 3. The modes belonging to the valley-diagonal irreducible 
representations correspond to the acoustic and optical modes at the original T point, according to which 
atoms of different sublattices oscillate in-phase or out-of-phase respectively. The modes belonging to the 
valley off-diagonal irreducible representations correspond to the modes at the I\± points, in particular, real 
linear combinations of displacements at both valleys. 
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Table 2: Classification of in-plane phonon modes according to the irreducible representations of Cq v . 
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Table 3: Classification of flexural phonon modes according to the irreducible representations of Cq v . 
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2.2. Theory of elasticity: graphene as a crystalline membrane 

The dynamics of low energy vibration modes of graphene are governed by its elastic constants, and 
a long-wavelength description in terms of a theory in the continuum is valid since we focus on in-phase 
displacements only. The mechanics of solids, regarded as continuous media, is the subject of the theory of 
elasticity [21, 22], which can be constructed by purely geometrical means with the only constraint imposed 
by the symmetry of the underlying crystalline structure. The unit cells of the graphene crystal can be seen 
as points on a surface embedded in R 3 , each of them labelled by a 3D vector r (throughout the text we 
will denote 3D vectors with upright bold font and 2D vectors with an arrow). Therefore, the elastic energy 
associated with such an object can be split into two terms with a different geometrical and physical origin, 
which we write as 

^elastic = Fat + -f\>) (2) 

where F st is a stretching energy term that accounts for the energy cost due to relative changes in the distances 
and in-plane bond angles between atoms, and iq, is a bending energy term which penalizes deviations from 
a flat configuration. 

The stretching energy depends on the first fundamental form, which is nothing but the pull-backed metric 
of the embedding space projected onto the surface, 

9ij = diV • djT, (3) 

where di, i = 1,2, is the derivative with respect to the internal coordinates which parametrize the surface. 
Note that this is a completely intrinsic quantity. In the Monge representation, carbon atoms at equilibrium, 
flat configurations are assumed to occupy positions labelled by ro = (cco, 2/o, 0) , where xq , yo correspond 
to the intrinsic coordinates that parametrize the surface. We consider deviations from these equilibrium 
positions given by a vector of displacements r = r 0 + u. The displacements are fields which depend on the 
position of the unit cell, u(xo,yo)- The first fundamental form is given by 

fjij = hij T 2 Uij , (4) 

where is the strain tensor. To the lowest order in the displacements, u = (u x . u y , h) , this reads 

u ij = ^ (diUj + djUi + dihdjh ) , i,j = x, y. ( 5 ) 

Note that Uij is quadratic to the lowest order in the out-of-plane displacements, h , as a manifestation of ah 
symmetry. We write the linear strain tensor only due to in-plane deformations as rt(ij) = ( d^Uj + djUi ) /2 
and the in-plane displacement as u = (u x ,u y ). The energy cost due to stretching can be written as a 
quadratic form in the strain tensor components [22, 21], 

F s t = \J d 2 rC^ kl u ljUk i, (6) 

where repeated indices are summed over and C is the elastic constants tensor. Since = Uji by 
definition, in principle there are only 6 independent elastic constants in a 2D solid. Additionally, the elastic 
energy density must be invariant under the point group symmetries of the crystal. Note that Eq. (6) is 
not a complete invariant under coordinate transformations since for that the integration measure should be 
completed with the factor y/detg « 1 + ua. The description in the continuum of a 2D hexagonal crystal is 
essentially different from a fluid or any continuous media in the sense that the action must not be invariant 
under diffeomorphisms given that the atoms in the solid define a special coordinate system [23]. 

For isotropic media, only two elastic constants are independent. Since in that case the response of the 
solid must be independent of the direction, the only tensor available to construct higher order tensors is the 
Kronecker delta <W. The only 4-rank tensors formed by satisfying the symmetries of are <Wc 5 kl and 
+ 5 ll 5^ k , therefore 

C ijkl = A 6 ij 6 kl + p (S ik 5 jl + 5 il S jk ), (7) 
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where A, /i are the Lame coefficients. In 2D. this is precisely the case of hexagonal crystals. In the particular 
case of graphene, the 3 independent components of u t j transform according to 


Itxx “b 'Uy y ^ A \5 

( ) ~ E - < 8 > 

Then, among the six independent combinations of the form UijUu, only two of them form scalars belonging 
to A\ irreducible representations. This is inferred from the algebraical reductions 

A\ x A± ~ Ai, 

A\ x E 2 ~ E 2 , 

E 2 x E 2 ~ Ai + A 2 + E 2 . (9) 

Therefore, for crystals with Cq v there are only 2 independent elastic constants corresponding to the Lame 
coefficients. Thus, the stretching energy for graphene has the explicit form 

Fat = ^ J d 2 f(Xuii + 2 nUijUij) . ( 10 ) 

Notice that inclusion of the quadratic term dihdjh in u.jj makes the theory defined by the above equation an 
interacting theory. For free membranes this has important consequences as will be seen in Section 3. If we 
discard these anharmonic terms in Eq. (10), we obtain two in-plane acoustic phonons with linear dispersion 
relations. At low momentum these modes can be classified as longitudinal and transverse, and the respective 
dispersion relations are given by 


= J^^\q\ = Vl \q\, ( 11 ) 

^ = \Zf ^ = vt ^ 

where p ~ 7.6 x 10 -7 kg/m 2 is the graphene’s mass density and Vl/t is the longitudinal/transverse sound 
velocity. Typical values for the elastic constants of graphene are /r ss 3A ss 9 eV A -2 [24, 25, 26], leading to 
v l/t ~ 10 4 m/s. 

Information about the bending of the membrane, or in other words, its embedding in 3 dimensional 
space, such as changes from point to point of the normal vector to the surface, n = d±r x d 2 ?/\dir x c? 2 r|, 
are described by the second fundamental form, 

Fij = n • didjv. (13) 

To lowest order in the displacement fields, and assuming smooth out-of-plane displacements, the normal 
vector at each point of the membrane is given by n ~ (—d x h, —d y h, 1). In this approximation the second 
fundamental form is given by 


Fij ~ didjh. (14) 

Two scalars may be constructed from the second fundamental form: the mean extrinsic curvature, F. n . and 
the intrinsic or Gaussian curvature detJy but note that the latter is a purely topological term which does 
not affect the equations of motion. The bending energy is then a quadratic form of the former [4], 

Fh = ^J ~ f / d 2 r(y 2 hf . 
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(15) 




Here k ~ 1 eV [26] is the (bare) bending rigidity. The dispersion relation of out-of-plane or flexural acoustic 
phonon modes of a free membrane is then quadratic, 



However, a rotational symmetry breaking tends to linearize the dispersion relation, for instance, due to the 
presence of strain. On the other hand, these modes are strongly suppressed in supported samples, where 
part of the spectral weight is transferred to a hybrid state with the surface Rayleigh mode of the substrate 

In¬ 
putting together the two contributions to the elastic energy of the membrane we obtain 

^elastic = \ J ( v2/l ) 2 + + Z/J-UijUij'j . (17) 

This equation has the same functional form as the usual von Karman free energy for thin plates [21]. 
However, two differences must be pointed out. First, for a thin plate the bending rigidity, k, is related 
to the elastic constants A and /i and the thickness of the plate. The same is not true for atomically thin 
membranes, for which it is not possible to define a thickness. Therefore, k is an independent parameter of 
the model for a membrane. Secondly, Eq. (17) corresponds to the potential energy of a Hamiltonian, not 
a free energy. The difference becomes important when thermal fluctuations, discussed in Section 3.1.1, are 
taken into account. For k sufficiently large, as in a plate, anharmonic interactions are suppressed and the 
free energy is given by the same expression as the potential energy. If n is small, the case of a membrane, 
then anharmonic interactions become dominant and the free energy can have a very different form than the 
potential energy. 

2.3. Electron-phonon coupling 

The symmetry approach has been widely applied to the problem of strained graphene ( e.g. [28],[29], [30]). 
With the help of group analysis it is possible to construct a systematic derivative expansion of the low 
energy effective Hamiltonian of graphene in the presence of non-uniform elastic deformations. For the sake 
of simplicity we will consider only the case of spinless graphene. The spinful case and the role of strain 
in spin relaxation are discussed in Refs. [31, 32], We restrict ourselves to linear order in the electron 
momentum k (first derivatives), so that the effective Hamiltonian will be a function of the strain tensor Uij 
and its derivatives, and the electron fields ip and ip . Beyond the Dirac approximation the procedure can be 
straightforwardly extended to higher order in k. In this way, the systematic expansion is characterized by 
two integers, ( n q , nu), the order of the derivatives of the strain tensor and the electron fields respectively. 

Since we are interested in long wavelength deformations that are not able to couple different valleys, 
in what follows we will construct the possible interactions allowed by symmetry based on the little group 
approach [33]. As explained in Ref. [34] it is sufficient to consider terms invariant under C^ v , the little 
group of one of the inequivalent Dirac points, together with the combined operation TC 2 , where C 2 is a 
7 T—rotation around the z axis. Once the degrees of freedom are established and the symmetries of the system 
are determined, all we need to know is the way they transform under the relevant group operations and how 
these can be decomposed as a sum of irreducible representations. Using standard group character operation 
[33] it is possible to determine how many and which are the terms allowed by symmetry at a given order 
( nq,n.k ). As explained in Ref. [34], up to first order in derivatives, n q + n*, < 1, there are only six terms 
involving the strain tensor. These six terms corresponding to K + are listed in Tab. 4 together with their 
physical significance and the corresponding relative sign at the inequivalent Dirac Point I\-. 

We can now write the most general Hamiltonian in the presence of non-uniform strain, up to first 
derivatives of the strain tensor (5), and up to first order in the derivatives of the electron fields ip and 
ip^ . It is important to note that Uij is a function of two independent variables Ui and h. In order to keep 
the necessary number of independent coupling constants, we will introduce for each H term, function only 
of u^, an extra term H, function only of h; 
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Hi 

(n q ,n k ) 

Interaction term 
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K 

Hi 

Hi 

(0,0) 

(0,0) 

0 U X x ^ }j yy )H 

i^xx 'U'yy')&x 2'U rC y(J 
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Position-dependent electrostatic pseudopotential 
Dirac cone shift or U(l) pseudogauge field , A y l ) 

+ 

h 3 
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Dirac cone tilt 
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H a 
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(u>xx 'Uyy){p’ x fox H - ®yfoy) 

Isotropic position-dependent Fermi velocity 
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h 5 

( 0 , 1 ) 
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Anisotropic position-dependent Fermi velocity 

+ 

H 6 

(1,0) 


dyiv^xx 'Uyy) ‘^^x'U'xy 

<Jz 

Gap opening by non-uniform strain 

- 


Table 4: Effective low energy terms in the Hamiltonian for the electron-strain interactions allowed by symmetry. The role of 
the gamma matrices as projectors into irreducible representations is clarified in [35]. 


6 6 

H = H Q + Y,9iHi + Y,9iHi. (18) 

i—l i=1 

The first term is the standard Dirac contribution Hq = Vp-(ia x k x + a y k y ), and the terms Hi are given in 
Tab. 4, and described below. The terms Hi are obtained from those in Tab. 4 through the substitution 
Uij —> dihdjh. Of course, the symmetry machinery can give us all the possible terms up to a given order in 
derivatives, but the coupling constants remain undetermined and have to be fixed with a concrete microscopic 
model. 

The second-quantized Hamiltonian operator is given by H = f d 2 rip' Hip, where the symmetric conven- 

J ^_ y 

tion for the derivatives acting on the electron fields is assumed, i.e., ip^kiip —> — i/2(ip^ di ip) = —i/2(ip^diip — diip^ip). 
In more general terms, dj acts only on the electron fields. The construction of second quantized hermitian 
operators becomes simpler using the symmetric derivative convention. Note that Tab. 4 gives the orders of 
the derivatives when terms are written with the symmetric convention. 

• Hi = ( u xx + u yy )l: This term has been already described in Ref. [36]. It is a scalar potential 

V s (f) = (• u xx (r) + u yy (r )). (19) 

The coupling g\ is called deformation potential and its strength has been estimated to be of the order 4 
eV for single layer graphene [37], although other works [38, 39] argue that screening leads to a strongly 
suppressed g\ ss 0. Its physical consequences have been explored in Ref. [40] and will be discussed in 
Section 4. 

• Hi = (u X x ~ u yy )u x — 2u xy a y : This term represents a shift in momentum space of the Dirac cone. It 
corresponds to the well known £7(1) elastic pseudo-gauge field 

A el {r) = ±52 {u xx (f) - Uyy (r ), -2 u xy (r)). (20) 

In the tight binding formalism 52 ~ P = — 9(logt)/i9(logo), where t is the hopping integral between 
nearest neighbours and a is the lattice constant. It is at the basis of most of the results in the literature 
related to strain engineering and there are experimental realizations [8] . It has also been used to explain 
data in artificial graphene [41], 

• H 3 = [{u xx — u yy )k x — 2 u xy ky] I: Dirac cone tilt. This term appears naturally in the description of 
the two dimensional organic superconductors [42] which are described by anisotropic Dirac fermions. 

It also arises when applying uniaxial strain in the zigzag direction, a situation that has been discussed 
at length in the literature [43, 44, 45, 37]. A novel physical manifestation of the Dirac cone tilt has 
been described recently in [46]. 

• H 4 = ( u xx + u yy )(a x kx + a y k y ): Isotropic position-dependent Fermi velocity [47]. 
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• H§ = UijtJikj ; i,j = x,y: Anisotropic position-dependent Fermi velocity [47]. This term, together 
with H 4 , has been predicted to arise within the geometric modelling of graphene based on techniques of 
quantum field theory in curved space [48] . It was later obtained in a tight binding model by expanding 
the low energy Hamiltonian to linear order in k and Uij [47, 49]. Since the Fermi velocity is arguably 
the most important parameter in graphene physics, this term affects all the experiments and will 
induce extra spatial anisotropies in strained samples near the Dirac point [50, 51, 52, 53, 54, 55]. This 
term will be described in detail in Section 2.4. 

• H e = [d y (u xx — u yy ) + 2 d x U xy ]<t ,: This is a very interesting term that suggests a new gap-opening 
mechanism [34]. It can be seen as the Zeeman coupling of pseudospin to the associated pseudomagnetic 
field B z = d x Ay l — d y A x [56]. The magnitude of this gap has been estimated to be 7 meV [34] . 

It is important to notice that because of the symmetric convention, some non-spurious terms seem to be 
missing, as it is the case of the vector field F, [47]. This is just an artefact of the symmetric convention and 
all the relevant terms will appear in the equation of motion. In the case of Tj, its presence becomes evident 
by integrating by parts and H$. 

2.3.1. Experimental confirmation of the elastic pseudomagnetic fields 

The term H 2 representing fictitious magnetic fields coupling with opposite sign to the two valleys was 
the first to be identified in the early publications and its physical consequences were discussed at length in 
the review [8]. The theoretical suggestion that Landau levels can form associated to strain in graphene was 
done in [57, 58]. The simple observation that applying a trigonal strain gives rise to uniform pseudomagnetic 
fields in a region of the sample, gave birth to real strain engineering. The theoretical proposal was soon 
followed by the experimental confirmation. A Landau level structure was observed in a STM experiment 
on graphene samples with triangular nanobubbles in [59]. Fields of up to 300 Tesla were estimated in the 
central region of the bubbles. This pioneer experiment was followed by many reports of elastic Landau levels 
in graphene samples [60, 61, 62], cold atoms [63], and artificial graphene [64]. Elastic gauge fields have also 
been described in bilayer graphene [65], other 2D semiconductors [66] and lately in the recently discovered 
three dimensional Weyl semimetals [67]. The effects of the elastic pseudomagnetic fields on various transport 
coefficients have been analysed in numerous publications [68, 69, 70] and will be described in more detail in 
Section 4. 

2.3.2. Experiments probing local variations in vf 

There are several ways to measure the Fermi velocity experimentally. In this section, we concentrate 
on measurements that are local and can thus probe the spatial variation of the Fermi velocity predicted by 
Eq. (26). The optimal experimental technique to do this is scanning tunnelling spectroscopy (STS), which 
is sensitive to the local density of states (LDOS) 1 . Since the LDOS depends on vf as p(E) = 2\E\/nVp, if 
the strain variation is smooth one may replace vf —t vf{x) and fit the slope of the measured LDOS to this 
expression to obtain a local measurement of vf- 

Two recent experiments have reported evidence of a spatially varying vf using this method. The first 
was performed on graphene grown on a Rh foil [71], which develops quasiperiodic ripples. Local variations 
of vf on different parts of the ripples were reported with a rather wide distribution of values, consistent 
with the strain induced by the ripples. In the second, graphene was grown on BN [72], where strained ridges 
were randomly formed. In a more quantitative comparison, the LDOS measured in these ridges reveals 
a spatially dependent vf which is shown to be correlated with the height of the ridges. This represents 
convincing evidence that vp fluctuations indeed originate from strain. 

An alternative way to measure vp locally is with Landau Level spectroscopy. In the presence of a 
magnetic field, the LDOS is rearranged into a series of well defined peaks at the Landau Level energies, 
given by E n = sgn(n)vF\/2ehBn. A measurement of vf can thus be obtained by fitting this expression as 
well. This type of experiment has been performed with randomly strained graphene as grown on Si02, and 
the measurement indeed reveals a local variation of vf of 5-10 % [73]. 


1 See section 5 for a discussion of how superlattices may also induce local variations of the LDOS that is seen by STM. 
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2.f. Quantum field theory in curved spaces 

One of the features that makes graphene unique is the way the electronic degrees of freedom couple to 
structural deformations of the lattice, and how this allows to modify its electronic properties in interesting 
ways. In the previous section, we have shown how the deformation of the lattice is described in terms of a 
strain tensor Uij , and what are the possible strain couplings allowed by the discrete symmetry of the lattice 
in the low energy Dirac Hamiltonian. In this section, we discuss a complementary approach to the study 
of the effects of deformations which provides insightful connections to the geometry of the problem. This 
approach originates in the idea that, since local distances in the strained lattice can be described in terms of 
a metric (the first fundamental form introduced in (3)), a natural model to describe the coupling of electrons 
to the deformation should be the Dirac equation in curved space. The formalism in full glory is taken from 
quantum field theory (QFT) in curved space [5] and was explained in detail in the review article [8]. In this 
section we describe the modern uses of the metric (unrelated to strain deformations) to generate response 
functions in condensed matter theory. 

While this approach is phenomenological rather than microscopic in nature, it offers an interesting 
geometrical interpretation of the different couplings, and it is a well understood formalism. The Dirac 
equation in curved space is completely determined by the metric gij = Sij + 2 mj and takes the form 

H = -i J <Fxy/g'ip'fea{di + tii)‘ l P, (21) 

where a = 1, 2 and * = 1,2 run over the space dimensions, */> = 7 ° 7 * = a 1 are the Pauli matrices, e\ are 

the tetrads, f1,; is the spin connection and ^Jfj is the volume factor. These three geometric quantities can be 
obtained from g- r] , see Ref. [47] for explicit expressions. Many authors have studied this model in connection 
with different problems, which include the spectrum of fullerenes [74], the modelling of disclinations[75, 76], 
graphene wormholes [77] proposals of metrics to generate in optical lattices [78], the holography conjecture 
[79] and the Unruh effect [80, 81, 82]. This model has also been applied to the study of the optical conduc¬ 
tivity in graphene [83, 84], to the scattering by smooth deformations in topological insulator surfaces[85, 86], 
and to study the Quantum Hall effect [87]. 

To be able to compare this approach with the standard tight-binding results, one proceeds by expanding 
Eq. 21 to first order in Uij for small strains, which can be written as 

H = -i[vij(x)<Tidj + o-jTi], (22) 

where 

Cjj = vp{6ij T SijUkk v.ij ) (23) 

Tj =2^3 v ij (24) 

where vf is the Dirac fermion velocity in the absence of distortion. The most interesting feature of this 
model is the presence of the matrix v. t j , which may be interpreted as a tensorial, spatially dependent Fermi 
velocity. This interpretation was first put forward in Ref. [48] , where it was shown that it can be responsible 
for local modulations of the density of states (DOS) that are correlated with strain. The model also displays 
a spin connection, which can be responsible for pseudo-spin precession, for example. It can be seen that 
the two terms in the expression for the Fermi velocity are consistent with those allowed by the discrete 
symmetry of the lattice, terms H 4 and H$ in Section 2.3. 

The model just presented has many nice features and has proven useful, but it has the drawback that it 
has not been derived from a microscopic Hamiltonian. The simplest microscopic derivation of the coupling 
of electrons to strain can be done with a tight-binding model, assuming that the hopping parameters change 
locally with distance, t n = t( 1 — g 2 \Au n \), where <72 is the dimensionless electron phonon coupling defined 
in the previous section. This model has been used extensively across the graphene literature to show that, 
to lowest order, strain couples as a gauge field of the form H 2 in the low energy Hamiltonian. The space- 
dependent Fermi velocity is also present in this model, when the expansion is taken to first order both in 
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strain and in momentum [47], which gives a Hamiltonian 

Htb = -i[vfj B (xfadj + v F cTiTi + iv F <JiAf], (25) 

where the field Af is the one given in the Hamiltonian H2 of Section 2.3, and is the Fermi velocity 
tensor obtained from tight-binding, which is 

vJ B = v F [Sij - + 5ijUkk)\, (26) 

with v F = 3/2to. Again, the two terms found from the microscopic approach are consistent with those 
found in the symmetry analysis, H 4 and H 5 . One may observe, however, that the coefficients obtained for 
these terms in tight-binding are not the same as those in the curved space model. The physical predictions 
of these two models are however qualitatively similar. The fact that the tight binding model gives the 
functional form of the Fermi velocity tensor that is different from the one postulated by curved space model 
was recently discussed in Ref. [88]. 

Analytical expressions for the dependence of the Fermi velocity on strain from which the coefficients of 
the terms H 4 and H 5 can be extracted were also reported in Refs. [89, 54]. The anisotropy in the Fermi 
velocity was demonstrated numerically in tight-binding studies [90, 91] and in ab initio calculations [37]. 
The corrections to the Fermi velocity tensor to second order in strain were computed in Ref. [92]. The 
approach of Refs. [93, 94] based on discrete differential geometry can also be used to compute higher order 
corrections. The tight binding model with space-dependent velocity has also been used to study bound 
states in strain superlattices [95] and to study the modulation of persistent currents in graphene rings[96]. 
The possibility of tuning the strain to get a given velocity profile has been suggested in [97]. 

It is worth emphasizing that the Fermi velocity tensor, as well as the other terms obtained from the 
symmetry analysis, was computed by expanding the Hamiltonian around the high symmetry point K. Only 
in this way one can classify the different terms under the symmetry of the little group at K. It has recently 
been argued that, in the case of constant strain, since the Dirac point is shifted due to the constant vector 
potential, the Fermi velocity tensor has a different form when the expansion is performed around the shifted 
Dirac point [93, 94]. 

There are several works that have addressed the modelling of a spatially dependent Fermi velocity with 
a different phenomenological model, which amounts to the direct replacement of v F by a scalar function 
v(x) in the Hamiltonian [98, 99, 100, 101, 102]. This model is equivalent to considering just the isotropic 
coupling to the strain H4 with a scalar function v(x) = v F Ukk- 

The motivation to study the case when the only term affecting the Fermi velocity is the isotropic one is to 
analyse how to control wave propagation by velocity changes in analogy with the modulation of the refraction 
index in optics. These works studied the propagation of Dirac electrons in different scalar velocity profiles, 
concluding that the transmission through a velocity barrier has a strong dependence on the incidence angle, 
but is always equal to one at normal incidence. When v F in the barrier is larger than outside it, there is a 
critical angle above which no transmission is possible, in analogy with the Brewster angle. 

2.4-1. Uses of the metric as an auxiliary field to generate response coefficients 

We here review some uses of the metric and the spin connection as auxiliary tools to derive exotic 
transport properties of materials. 

In an effective action framework transport coefficients can be computed from the partition function of 
a system by coupling appropriate source terms to the operators whose response we want to analyse. The 
best example is that of the electromagnetic currents obtained by coupling a background electromagnetic 
potential (even in the absence of actual fields). Then Kubo formulae can be used to compute the coefficients 
of the current-current correlation functions. 

Thermal transport is one of the most difficult aspects to address in materials. In a seminal work J. M. 
Luttinger [103] proposed the use of the metric as an auxiliary tool to analyse thermoelectric transport co¬ 
efficients in solids. In the remarkable footnote 7 he writes: In fact, if the gravitational field didn’t exist, one 
could invent one for the purpose of this paper. The aim was to derive an expression for the thermal conduc¬ 
tivity similar to the Kubo formula of the electrical conductivity. The main difficulty with thermal transport 
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coefficients is that there is no Hamiltonian which describes a thermal gradient since the temperature is a 
statistical property of the system. Inspired by the relation between the metric and the energy-momentum or 
stress tensor with Greek indices running over space-time coordinates) that holds in general relativity: 

T^ = ^=^ S -. (27) 

V~9 d 9»» 

(S is the action of the matter) he proposed to couple the electron system to an auxiliary field (an inho¬ 
mogeneous gravitational field) which causes energy or heat currents to flow. The interaction term was 
Sint = / h(r)$(r,t) where h{r) is the Hamiltonian density of the matter system and <h(r, t) is a “gravita¬ 
tional potential”. A varying $ will give rise to a varying energy density which in turn will correspond to a 
varying temperature. No matter how ambiguous the original work can look today, it paved the way to very 
important recent developments in condensed matter. An attempt to put this work on firmer grounds by de¬ 
riving expressions similar to (27) for a non-relativistic system has appeared recently in [104]. The proposed 
formalism involves the Riemann Cartan geometry, which, in addition to curvature, includes torsion in the 
background geometry (for a "primer” on torsion in condensed matter applications see [105]). The viscosity, 
another important transport coefficient of the electronic fluid, also involves the metric in its definition as 
will be discussed in the next subsection. 

2.5. Topological aspects and time-dependent deformations 

Topological phases of matter are a new paradigm in condensed matter physics. These phases add to the 
standard classification in terms of local order parameters and broken symmetries new quantum numbers 
associated to topological invariants [106,107] . Since the early developments of the field and due to its peculiar 
properties, the honeycomb lattice has been used as the ideal model to demonstrate topological properties. 
The pioneer works proposing this lattice as a “condensed matter simulation of a three-dimensional anomaly” 
[108, 109], together with the analysis of the spin orbit in graphene done in Refs. [110, 111] opened the modern 
field of topological insulators. Most of the new topological materials (graphene, topological insulators and 
superconductors, Weyl semimetals [112, 113]) share the property that their low energy electronic properties 
are described by Dirac fermions [114, 115] in one, two or three spacial dimensions enlarging and widening 
the bridge between high energy and condensed matter in this century. 

At a zeroth order approach to topological aspects, Dirac points in (2+1) dimensions are usually protected 
by a winding number: the Berry phase acquired by the spinor wave function when circling around a Dirac 
point in momentum space. Time reversal symmetry implies that Dirac points in crystals arise in pairs with 
opposite winding numbers [116]. Lattice deformations move the Dirac points that eventually can merge 
giving rise to interesting topological phases. The existence and merging of Dirac points in crystals was 
studied in the early work [117] and a very complete analysis of the issue has appeared in the recent review 

[115]- 

In this section we review the late developments on the influence of lattice deformations on the topological 
properties of graphene and related materials. The topological phases induced by a combination of strain 
and interaction will be discussed in Section 4.5.1. 

2.5.1. Hall viscosity 

In a similar way as done by Luttinger with the heat transport, the metric is being used more recently 
to analyse properties such as the Hall viscosity [118] in topological insulators and in quantum Hall fluids 
[119, 120, 121, 122]. The analogue of friction in fluid dynamics is provided by viscosity, which causes 
dissipation of the energy (momentum) flow. In elasticity theory we can define (see Section 2.2) 

T ij = C ijkl u kt + r] ijkl u k i, (28) 

where Uij is the strain tensor, and C l ^ kl and are the elastic and viscosity coefficients. Viscosity can 
then be defined as the response of the system to a time-varying strain. In general, this viscosity tensor 
possesses both symmetric and antisymmetric components under the permutation of pairs of indices. While 
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the symmetric part is generally associated to dissipation and vanishes at zero temperature, the antisymmetric 
part arises when time reversal symmetry is broken and was first described in the quantum Hall fluid [118]. 
A stress-stress form for the response function yields a Kubo formula for the viscosity. 

The Hall viscosity can be derived by coupling the system to an external metric gij. Under a small 
time-dependent shape deformation which preserves the volume: gij = Sij + Sgij(t ), the stress tensor is 

T ij = PS ij - \n ijkl d t 6 gij , (29) 

where P is the pressure and the tensor ijijki = —gkiij is the Hall viscosity. This transport coefficient is 
particularly interesting in the case of fractional quantum Hall fluids (FQHF). Unlike the Hall conductivity 
which is quantized and dimensionless, the Hall viscosity has units of l/L 2 . In a quantum Hall system the 
Hall viscosity, defined as the response of the electronic fluid to geometric shear distortions, is proportional 
to the density of intrinsic angular momentum: tja = 1/2 sn (n is the particle number density and s is the 
average spin per particle, so sn is the area density of intrinsic angular momentum). Computation of the Hall 
viscosity in topological matter and its relation to the total angular momentum is one of the most active areas 
of research nowadays involving all the geometrical tools concerning the study of QFT in curved background 
fields including torsion [119, 122, 123]. In a very recent publication [124], Ward identities have been used to 
relate Hall viscosities, Hall conductivities and the angular momentum. 

In the early publication [125] it was recognized that Hall fluids couple also to the curvature of the 
background surface, a fact that manifests itself by the presence of a shift in the relation between the number 
of electrons N e , and the number of magnetic flux quanta AU going through the manifold. In a flat system 
they are related by = N e . where v is the filling factor. In a compact curved manifold such as a sphere 
there is a shift S such that: 

N.,, = v~ l N e - S. (30) 

It was observed that the shift on a sphere and in a torus were different (it was zero in a torus) what 
led the authors to conjecture that the Hall fluid did couple to the spin connection vector form defined as 
= uj^Eab- Coupling this vector as a source in an effective action formalism allows to obtain the shift as 
the coefficient of a mixed Chern-Simons electromagnetic-gravitational term of the form Swz ~ £ l ^ k uJ i 'S/jOk, 
known as the Wen-Zee term. The effective Chern Simons action describing the geometrical factors in Hall 
fluids has been put in firmer grounds recently in [123]. 

The Wen-Zee action can also be seen as an example of anomaly related transport, a field that is also 
investigated very actively these days in both condensed matter and high energies in the context of the quark- 
gluon plasma [126]. The most commonly cited example of the new non-dissipative transport phenomena 
occurring in the quark-gluon plasma is the chiral magnetic effect [127] that refers to the generation of an 
electric current parallel to a magnetic field whenever an imbalance between the number of right and left- 
handed fermions is present. Another interesting example is the axial magnetic effect (AME) associated with 
the generation of an energy current parallel to an axial magnetic field, i. e. a magnetic field coupling with 
opposite signs to right and left handed fermions. These new transport phenomena are specially relevant 
to the (3+1) dimensional Dirac materials, particularly in Weyl semimetals. The axial magnetic effect in 
Weyl semimetals as a fingerprint of the gravitational anomaly has been proposed in [128]. Of particular 
interest is the part associated to the anomalous transport coefficients generated by the gravitational anomaly 
[129, 130, 131], 

2.5.2. Transverse conductivity from time-dependent deformations 

The valley degree of freedom in graphene has been proposed as a novel ingredient to encode information, 
its control and manipulation led to the concept of valleytronics [132]. As mentioned in Section 2.3, lattice 
deformations and defects couple to the electronic current in graphene in the form of vector fields that couple 
with opposite signs in the two valleys. This fact has given rise to several interesting proposals related to 
valley physics in strained graphene [133]. It was shown in [50] that the combination of real and fictitious 
magnetic field induces a valley asymmetry in the Landau level structure. 
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Figure 2: (colour online) Comparison between the Hall effect (upper), the valley Hall effect (middle), and the effect proposed 
in the text (lower) following the scheme of the original implementation of the experiment first proposed by Hall. The arrows 
indicate the flow of valley polarized electrons under the action of an external electric field or voltage (upper and middle) and of 
a time dependent strain (lower part). In the later case the total current is < J X K + J X K , E? 1 .(Fig. adapted from ref. [135]). 

Following Eq. (20), a time-dependent deformation gives rise to a “synthetic electric field” Ef = <9 0 Af 
that also couples with opposite signs to the two Fermi points. The observable consequences of the “synthetic 
electric field” in graphene were described in the early publication [134]. It was recognized that, due to the 
opposite signs of the coupling at the two valleys, this synthetic electric field can drive charge-neutral valley 
currents in carbon nanotubes and graphene that are not subjected to screening and hence can have a stronger 
influence on the electronics than the deformation potential. 

An interesting topological response due to the interplay of electromagnetic and a time-dependent elastic 
vector field was studied in [135]. The total Chern number - and hence the Hall conductivity - of a time- 
reversal invariant system is zero. The topological signals in these systems, such as the spin Hall effect, have 
associated dissipationless transverse currents involving degrees of freedom - spin, valley, layer [136] other 
than the electric charge. It was shown in [135] that the interplay of electromagnetic and elastic vector fields 
can give rise to a transverse electric current. The argument goes as follows: In the presence of an external 
electromagnetic and elastic field the gapped graphene Hamiltonian reads: 

K = Y^^l,%( T(T * k *+ a V k v)'l > T,% ~ Yy^r,k T(Jx ( eA x m + Tg 2 Al l ) V- T) fe] - 

r r 

- Mr + T 92A e y l ) a^ rX ] , (31) 

T 

where r = ±1 denotes the valley degree of freedom, e is the electric charge, A em and A el stand for the 
electromagnetic and the elastic vector fields respectively, and the parameter g 2 > 0 encodes the strength of 
the elastic coupling. The coefficient r in front of the coupling of the electronic current to the elastic vector 
field marks the difference with the electromagnetic gauge field that couples with the same sign to the two 
valleys. In the presence of a time-reversal preserving mass m (having the same sign in both valleys) as the 
one considered in (31) there is a non zero Chern number which takes opposite values at the two Fermi points 
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Cg = —Cg, = sgn(m)/2. In the absence of elastic deformations the response to an external electric field 
E em is a transverse charge current with opposite signs at each Fermi point: (J*) = e 2 C T e l i Ef m , so the total 
charge current (Jb + Jb ( ) vanishes. The synthetic elastic field will induce a charge response in the system 
at each Fermi point: (J 2 ) ~ rC T e' < ■ } E| , . Now, because Cg = —Cg, the total net charge current is non zero 
and its value is twice larger than the value at each Fermi point: 

(J i ~+ ji i}l )~2C R e ij Ef. (32) 

This result can also be understood from the mixed Chern-Simons action arising from the combined vector 
field Vi = Af fn + rAf . Within a functional integral approach one can integrate out the fermionic degrees of 
freedom f> T g, ijy - in the action derived from Eq. (31) and write the odd part of the effective Lagrangian: 

C eff = 2 eg 2 Cge^A; m d p At l + 2 eg^CAp p A™ + 

+ + J£A; 1 , (33) 

where we have added to the Chern Simons action the external sources: J M is the total charge density current 
that naturally couples to the electromagnetic field, and ,/(.) is a classically conserved current associated to 
the elastic held AT},. Notice that the standard Chern-Simons term bilinear in A em or A el vanishes for the 
total action due to the opposite value of Ck at the two valleys. Only the mixed term survives. The total 
charge current density ( S e jj = f d 3 C e g) is: 

P‘> = Mr = 2 ^ c « e ” A f s 2eg ^f' lE f- (34) 

where for simplicity it has been assumed that A} = 0 and have replaced £ lj0 = . Notice that the 

generation of this transverse electric current in a time reversal invariant system does not contradict the 
general result of ”no Hall current in T-invariant systems” because the Hall current is the response to an 
electric Held while the present current is the response of the system to a time-dependent elastic deformation. 
A schematic description of the various possible valley currents is presented in Fig. 2. The action discussed 
here is similar to the Wen-Zee action described in Section 2.5.1 but there the additional vector Held that 
couples to the electromagnetic field in the mixed Chern Simons term was induced by the geometric spin 
connection (higher order in derivatives, see discussion in Section 2.3) while here it is the elastic gauge field. 
The result in both cases is similar: an electromagnetic response induced by a time-dependent deformation. 
This can also be seen as a quantum analogue of standard piezoelectricity. 

Elastic electric fields have also been analysed recently in [137, 138, 139]. As it is known, what prevents 
graphene from having observable topological properties is the the small value of the spin orbit coupling and 
the difficulty to open a gap in the spectrum. It has recently observed that depositing graphene on hexagonal 
BN opens a sizeable gap and a valley Chern invariant arises in commensurate stacking [140]. These type of 
heterostructures are discussed at length in Section 5. 

3. Lattice anharmonic effects in crystalline membranes 

3.1. Theory of anharmonic crystalline membranes 

3.1.1. Classical theory for ideal free, flat membrane 

The usual starting point for the study of lattice properties of a solid is the harmonic theory for lattice 
vibrations. The harmonic theory is valid provided that the fluctuations of the atoms around the equilibrium 
positions are much smaller than the average interatomic distance. In bulk three-dimensional crystals, the 
harmonic theory, or minimal extensions of it as the quasi-harmonic approximation, is usually a sufficient 
approximation for most applications. However, the situation is quite different in two-dimensional crystalline 
membrane (membranes with a finite shear modulus, also referred to as solid membranes), such as graphene 
and other two dimensional crystals, where, due to the low dimensionality, atomic fluctuations are not small. 
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As a matter of fact, it was argued by Peierls [141] and Landau [142] that two dimensional crystals should not 
exist at any finite temperature. Mernrin [143] formalized this argument and proved that, provided the pair 
interaction between two atoms decays faster than r -2 (r is the distance between atoms), there is no positional 
long-range order. However, there is a residual quasi-long range order, characterized by correlations with a 
logarithmic behaviour with distance, which is typical of two dimensional systems. Nevertheless, Mermin 
also pointed out that orientational long-range order can still exist. 

For a free two dimensional crystalline membrane embedded in a three dimensional space, the situation 
is potentially worse since now the membrane has the possibility of fluctuating in the third dimension. 
Therefore the study of properties of the flat phase in crystalline membranes, its thermodynamic stability 
and the study of in-plane and out-of-plane positional and orientational orders became a topic of intensive 
research in the context of statistical mechanics (for an in depth discussion on the physics of membranes, the 
interested reader can refer to Refs. [4, 144, 145]). It turns out that, although a free crystalline membrane does 
not have positional long range order, the anharmonic interaction between in-plane and out-of-plane modes 
is essential in reducing fluctuations and protecting the long-range in-plane and out-of-plane orientational 
orders, such that a flat configuration of the membrane is thermodynamically stable. [3, 23, 146]. 

Much of the physics of flat crystalline membranes can be understood using the simple classical model 
of continuous local elasticity with a potential energy described by Eq. (17). [3, 23] As already discussed 
in Sec. 2.2, the inclusion of the term dihdjh in the strain tensor Uij makes the theory defined by Eq. (17) 
interacting. Neglecting anharmonic effects, in-plane and out-of-plane modes are decoupled and at classical 
level the theory is described by the correlation functions 
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where P^ T (q) is the longitudinal/transverse projector along the momentum q ., ujS, uk and ufk are, respec¬ 
tively, the dispersion relations for the flexural, in-plane longitudinal and in-plane transverse phonons (see 
Eqs. (11), (12) and (16) of Sec. 2.2), and ()o is a thermal average with respect to the harmonic theory. Using 
the harmonic correlation functions to compute the on-site displacement-displacement correlation functions 
one obtains 


(“ 2 )o = j 0 ~log(L/o), (37) 

(h 2 )o = J o ^ L 2 , (38) 

where L is the linear size of the membrane and a is an interatomic distance. The logarithmic dependence of 
(u 2 ) q with the system size is the expected behaviour of quasi-long range order in two dimensional systems. 
However, the linear growth of the out-of-plane fluctuations with system size, yj ( h 2 ) 0 ~ L, indicates that 
at harmonic level a flat membrane is not possible, and that the membrane is instead crumpled. There is 
however a long range in-plane orientational order, which can be seen by looking at [143] 


(3jUj Ok a i j q 



qiQk(uj(q)ui(-q)) 0 , 


(39) 


which remains finite. There is also a quasi-long range order of the local normals of the membrane, n(x), 

(<Sn-(5n> 0 ~ J ^J^\q\ 2 (h(q)h(-q)) 0 ~ log (L/a), 
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where, to lowest order in the displacements, Sh(x) = n(x) — z ~ — (V/i(x), 0), with z the unit vector 
pointing in the reference out-of-plane direction. Nelson and Peliti [3] were the first to argue that the 
interactions between the in-plane and out-of-plane fluctuations present in Eq. (17) would give origin to 
a long-range interaction between local Gaussian curvatures (which to lowest order in the displacements 
reads —S7 2 P^dihdjh/ 2), that would stabilize the flat phase. Since Eq. (17) is quadratic in the in-plane 
displacements, zq, these can be integrated out, and the effective action describing the dynamics of the 
out-of-plane mode is given by [3] 

S e g[h] = ^ J d 2 r(d 2 h) 2 + ^ J d 2 f(P^dihdjh) 2 , (41) 


where Y = 4/x(A + /z)/(A + 2/z) is the two dimensional Young modulus of the membrane. Computing to 
lowest order in perturbation theory the self-energy correction to the out-of-plane mode correlation function, 
{h(q)h(—q}} = ksT (n\q\ A + £(q)) \ due to the anharmonic terms in Eq. (41) one obtains [3, 4, 147, 148, 
149] 


< 42 > 

As we go to low momentum we have that T,(q) /t|< 7 ] 4 , which shows that perturbation theory breaks down 

and that a crystalline membrane is a strongly interacting system. Using a Ginzburg like criterion, from 
the condition £(g) ~ ft|g| 4 , the momentum scale below which anharmonic effects become dominant is given 
by [3, 4, 150, 147, 148, 149] 
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(43) 


which for graphene has the value q* ~ 0.2 A^ 1 . This momentum scale can be translated into a length 
scale L* = 2n/q* — 3 nm, that is typical size of a membrane above which anharmonic effects become 
dominant. For large membranes and below the anharmonic scale q *, Nelson and Peliti [3] showed that a 
non-perturbative treatment calculation of the self-energy leads to a correlation function for the out-of-plane 
displacement characterized by an anomalous momentum dependence. A renormalization group calculation 
[23] showed that the in-plane correlation functions also acquire an anomalous momentum dependence. In 
general, at small momenta the behaviour of the displacement correlation functions is characterized by two 
exponents rjh and r) u , such that 


(h(m-q))~Q-* +Vh , ( 44 ) 

(ui{q)Uj(-q)) ~ q^ 2 -^. (45) 

This behaviour can also be seen as the effective elastic constants becoming momentum dependent at low 
momenta 


K e s(q) ~ q Vh , (46) 

Aeff(g),Meff(g) ~ q r,u - (47) 

Using the renormalization group method, it was shown that the exponents r/h and rj u are not independent, 
but are related by [23, 146] 

Vu + 2Vh = 4 - D, (48) 

where D is the dimension of the membrane, with D = 2 for physical membranes. Provided r/h > 0 and 
r] u < 2 (for D = 2, these are equivalent conditions, provided Eq. (48) holds), there is in-plane and out-of¬ 
plane long-range orientational order, since taking into account the momentum dependence of the effective 
elastic constants yields 


(49) 

(50) 
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{diUjdkUi) ~ L Vu 2 , 
(<5n • Sh) ~ L~ Vh , 





and therefore a flat phase for the membrane is indeed stable. It was argued in Refs. [151, 146] that the lower 
critical dimension above which a flat phase for a crystalline membrane is possible should be smaller than 2. 
Anharmonic effects also reduce the crumpling of the membrane, since the height fluctuations scale with the 
system size as 

(h 2 ) ~ L 2 ^ (51) 

where £ = l—rjh/2 is the roughness exponent of the membrane. At the same time, for r] u > 0, the anharmonic 
effects generated by the out-of-plane fluctuations lead to a degradation of the positional quasi-long range 
order of the membrane, since 

(u 2 )~L r '“. (52) 

The anomalous dependence on momentum of the effective elastic constants was confirmed using the self- 
consistent screening approximation for the continuous model described by Eq. (17) [152, 153, 147], Monte 
Carlo and molecular dynamics simulations [154, 155, 156, 157, 158, 159] and functional renormalization group 
methods [160, 161, 162, 163]. The relation between anomalous exponents Eq. (48) was also derived using 
the self-consistent screening approximation [164, 152] and was further tested using Monte Carlo simulations 
[155, 156]. 

Although it is believed that the anomalous exponents should be universal quantities, there is some disper¬ 
sion on the reported values. An analytical result obtained using the self-consistent screening approximation 
(which becomes exact in the limit of large codimension of the membrane, dc , which for a physical membrane 
if given by dc = 3 — 2 = 1) gives a value for r] h of 0.821 [152], an improved self-consistent calculation to next 
leading order in 1 /dc gives a value of 0.789 [153], functional renormalization group calculations give 0.849 
[160, 163] and 0.85 [161, 162], molecular dynamics simulations give 0.81 [154], Monte Carlo simulations give 
0.72 [155, 156], 0.85 [157] and 0.795 [158], while renormalization group Monte Carlo gives 0.79 [159]. 

3.1.2. Deviations from the ideal case: defects, corrugations and residual strains 

The previous section concerned ideal membranes without defects and unconstrained by any external 
stress or fixed boundary condition. In an experimental setup graphene samples are either supported on a 
substrate or suspended over a trench. In supported samples coupling to the substrate leads to the opening 
of a gap in the out-of-plane phonon dispersion relation and, therefore, anharmonic effects should be severely 
quenched [27]. A suspended sample is much closer to the ideal case. Nevertheless, even a suspended sample is 
still subject to constrained boundary conditions, residual stresses, impurities and other sources of disorder. 
Therefore, it is important to understand how the previous discussion is affected once we consider these 
effects. 

The possibility of buckling due to constrained boundary conditions was studied in Refs. [165, 166]. The 
effect of residual stresses was studied in Ref. [167], where it was found that the anomalous momentum 
dependence of the correlation functions survived for small externally applied stresses. It was estimated that, 
in graphene, an applied stress corresponding to a strain of 0.25% can completely quench anharmonic effects. 

The effect of disorder in the stability of flat membranes was studied in Refs. [168, 169, 170, 171], where 
disorder was simulated by random strain and curvature fields. It was found that disorder can indeed 
destabilize the flat phase and that this instability can still be present at finite temperature, depending on 
the type of disorder. This possibility was verified with atomistic simulations of graphene [172], where it 
was found that a 20% adsorption of OH molecules leads to crumpling of the membrane. Membranes with 
a frozen non-flat background metric characterized by out-of-plane variances with a power law dependence 
on momentum, called warped membranes, were studied in Refs. [173, 174]. In these works it was found 
that for strong enough fluctuations of the background metric the anomalous exponents can differ from the 
ones due to thermal fluctuations in a flat membrane. The effect of short range quenched curvature disorder 
was also recently studied in Ref. [175] using a renormalization group analysis. There it was found that the 
experimental results of Refs. [176] and [177], for the average size of ripples and the mean deviation of the 
normals in suspended graphene samples, compare well with values obtained with the renormalization group 
analysis for weak curvature disorder. 
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3.1.3. Quantum theory 

The previous discussion was based on a classical description of membranes. The original theoretical works 
on crystalline membranes were motivated by the study of biological amphiphilic membranes [4] . Since such 
objects are formed by very large molecules, one can expect that their mechanical behaviour will be governed 
by classical statistical mechanics at room temperature. The same assumption is not easily justifiable in 
atomically thin membranes such as graphene. A simple estimation of graphene’s Debye temperature, gives 
us Td ~ 1000 K, considerably larger than room temperature. In this scenario one expects that quantum 
effects will still be important at room temperature. 

The quantum treatment of flexural phonons was first performed in Ref. [178], in the context of the 
study of the combined effects of phonon-phonon and electron-phonon interactions in graphene. Performing 
a self-consistent calculation, it was found that at zero temperature and in the absence of electron-phonon 
interactions, anharmonic effects lead to logarithmic corrections to the bending rigidity of graphene. The 
effect of quantum fluctuations at zero temperature on the Young modulus was studied in Ref. [179] using 
an expansion in the number of out-of-plane directions of the membrane (which for a physical membrane is 
one). It was found out that quantum fluctuations contribute to a reduction of the Young modulus at small 
momentum scales, just as thermal fluctuations, although being a much smaller effect. It was estimated that 
thermal fluctuations are dominant with respect to quantum fluctuations for momentum scales smaller than 
the thermal wavelength for flexural phonons 



Anharmonic crystalline membranes at zero temperature were also studied in Ref. [180] by performing a 
renormalization group study. It was found that zero temperature fluctuations lead to an increase of the 
effective coupling constant between in-plane and out-of-plane phonons which leads to a destabilization of 
the flat phase of the membrane. This behaviour is to be contrasted with the effect of thermal fluctuations, 
which lead to an increase of the bending rigidity at large length scales, and therefore stabilize the flat 
phase. Refs. [178, 179, 180] did not take into account retardation effects for the in-plane phonons, effectively 
treating them as classical, on the basis that typical frequencies of in-plane phonons are much larger than the 
frequencies of out-of-plane phonons at low momenta. The effect of retardation of the in-plane phonons was 
studied in Ref. [148, 149] by means of a perturbative calculation of the out-of-plane phonon self-energy. It 
was found that such retardation effects, change the corrections to the bending rigidity from being logarithmic 
to a power law in momentum, contributing to an increase of the bending rigidity at long wavelengths, a 
scenario that is very similar to the one due to thermal fluctuations. Furthermore, it was estimated that the 
effect of quantum fluctuations is dominant for temperatures smaller than 


2 h (A + 2n) 2 
~ VW Y 


(54) 


where / is a dimensionless cutoff dependent quantity. For graphene it is estimated that T* ~ 70 — 90 K. 
The momentum scale below which anharmonic effects dominate, in the low temperature regime, was found 
to be smaller but of the same order of magnitude as in the high temperature limit Eq. (43). Quantum 
corrections to thermodynamic and elastic properties of graphene were also studied in Ref. [181], by means of 
a self-consistent calculation of the density matrix, where it was found that quantum and anharmonic effects 
become comparable for temperatures smaller than 100 K, while quantum fluctuations still play a significant 
role for temperatures as high as 1000 K. 

The contradictory nature of the results published so far does not allow one to access what is the exact 
role and importance of quantum effects in the physics of crystalline membranes, such as graphene. 


3.1.4- Anharmonic effects and electron-phonon interaction 

As we have seen in the previous sections, the flat phase of a crystalline membrane is very fragile and 
anharmonic effects play a very important role. In membranes that have low energy electronic excitations, 
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such as graphene, one might ask if these electronic degrees of freedom can have a significant impact in the 
physics of crystalline membranes. 

This question was first addressed in Ref. [182]. There it was argued that the deformation potential 
interaction between graphene electrons and phonons, Eq. (19), gives origin to an effective strain of the form 
a = gin, where n is the electronic carrier density. This coupling is asymmetric with respect to electron or 
hole doping. While for electron doping, n > 0, the membrane would be subject to a tension, in the case 
of hole doping, n < 0, the membrane would be subject to a compressive strain, which would lead to the 
buckling of the membrane. Such a structural deformation induced by electron-phonon interaction bears 
some resemblance with the physics of the Peierls transition and formation of charge density waves. The 
role of electron-phonon interactions in undoped graphene was studied in Ref. [183], where it was shown 
that the decay of flexural phonons into electron-hole pairs, leads to the disappearance of the flexural mode 
branch below a critical momentum. It was speculated that such disappearance could be associated with the 
formation of ripples. A perturbative evaluation of the normal-normal correlation function in the classical 
limit performed in Ref. [53], also displayed a peak at finite momentum induced by the electron-phonon 
interaction. Such peak was interpreted as signalling the onset of ripples. The same problem was studied 
at zero temperature in Ref. [178], via a self-consistent evaluation of flexural phonon correlation function 
and it was found that for large values of the deformation potential, g\ ~ 23.1 eV, the effective bending 
rigidity vanishes at a finite momentum scale of the order of q ~ 0.1 A. Similar conclusions were found in 
Ref. [179], where, by treating electrons and phonons on equal footing, it was found that electron-phonon 
interaction can drive the system towards an instability in the charge carrier density - Gaussian curvature 
channel at finite momentum, that was speculated to be associated with the spontaneous and simultaneous 
formation of structural ripples and charge puddles. This mechanism is suppressed by the presence of an 
electronic bandgap, and should therefore be more relevant in graphene than in semiconductor transition 
metal dichalcogenides or boron nitride. The nature of this electron induced instability was further studied 
in Ref. [184], by means of a self consistent calculation and a renormalization group calculation. It was 
found that the vanishing of the effective bending rigidity and the condensation of the mean curvature 
(which is given by d 2 h) are simultaneous and manifestations of the same transition at finite momentum. 
This transition should occur even in the absence of an external tension and does not involve an in-plane 
distortion. Nevertheless, the nature of the relevant operators at this critical point is still not clear, and a 
condensation of the Gaussian curvature could not be ruled out. Therefore, further analysis is needed to 
clearly identify the electron-phonon interaction as a candidate to the formation of ripples in graphene. 

3.2. Mechanical properties 

3.2.1. Structure of suspended two dimensional crystals 

One of the most important quantities in the study of the structure of matter is the static structure 
factor S(q) = £ ( e lc G R "- R ™)), where R„ are the atomic positions (in 3D space), q = (q x , q y , q z ) is the 

transferred momentum and the average is to be understood as a time or ensemble average. As discussed 
in Sec. 3.1.1, in a crystalline membrane displacement-displacement correlation functions have an anomalous 
momentum dependence, which will affect the structure form factor. In a scattering experiment, Bragg peaks 
will no longer be sharp, but will instead be broadened with an algebraic decay. It is in this weaker sense, 
that a two dimensional crystal is defined [9, 185]. The structure factor of rough, self-affine membranes 
and surfaces 2 , characterized by power law displacement variances: (^\h(x) — h(x')\ 2 ^ = A \x — x'\ 2< ' and 

(\u(x) — u(x')\ 2 ^i = B\x — x'\ r, ' t , has been theoretically studied in several works [186,187, 164,188,189,190]. 
A semi-analytical expression for the structure factor of flat membranes subject to thermal fluctuations was 
obtained in Refs. [187, 164, 188]. For q = 0, (with q = ^ jq 2 + q^) it was found that the structure factor 

has a power law behaviour S(q = 0,q z ) ~ q z , while S(q,q z = 0) did not show a well defined power 


2 A surface is said to be self-affine if its roughness is described by a random variable with probability distribution that is 
invariant under anisotropic scaling transformations of the form x —> bx and z —)► b^z 
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law behaviour. In Ref. [189] the orientationally averaged structure factor, S(|q|), was studied in detail. It 
was found that for 1/L <C q <C l/£, where L is the membrane size and £ is the characteristic length scale 
of out-of-plane or in-plane fluctuations, one obtains S(|q|) ~|q|~ 2 as i n a completely flat, non-fluctuating 
membrane [186]. For l/£ «g<C 1/a, where a is the lattice spacing, it was obtained 5(|q|) ~|q| < ’ -3 , for 
the case where in-plane fluctuations are small compared to out-of-plane, and jS(|q|) ~|q| 2/C— 2 , j n the 
opposite case. 

The roughness exponent, (, of red blood cell membrane skeletons was obtained from the structure 
factor measured with light and x-ray scattering [191]. The roughness exponent was extracted by fitting the 
orientationally averaged structure factor to the power law behaviour S(|q|) ~|qp -3 , having been obtained a 
value of ( = 0.65 ± 0.10, which is in good agreement with the theoretical prediction for thermal fluctuations 
in a flat crystalline membrane (see Eq. (51) and text below). The exponent rjh was also measured in 
amphiphilic films [192], by fitting the structure factor measured with x-ray scattering, with a continuous 
model that neglects in-plane fluctuations. The best fit was obtained for r/h = 0.7 ± 0.2, once again in good 
agreement with the theory of thermal fluctuations in a flat membrane. 

It was experimentally found that suspended samples of both graphene and bilayer graphene display 
roughness [176, 193]. Transmission electron microscopy revealed that samples with lateral sizes of ~ 500 
nm display broadened Bragg peaks with Gaussian shape. These roughness showed no preferred orientation, 
and is associated with static ripples. It was estimated that these ripples have a characteristic lateral size 
of 2 - 20 nm and a mean deviation of the normals of the membrane of 5% in single-layer and 2% in bilayer 
graphene. The typical height of the ripples was estimated to be of 1 nm. Similar corrugations have also been 
observed in suspended samples of M 0 S 2 [194]. The roughness of suspended graphene was further studied 
using low-energy electron microscopy and diffraction in Ref. [195]. It was found that at small length scales, 
suspended graphene appears to be a self-affine surface characterized by a correlation length of £ ~ 24 nm 
and a roughness exponent of £ ~ 0.5. At larger length scales, graphene seemed to present a preferred ripple 
wavelength of A r i pp i e ~ 60 nm. Existence of the two length scales £ and A r j pp i e indicate that at large length 
scales graphene deviates from a purely self-affine surface and seems to be a mounded surface. [195] It was 
also found that the roughness exponent becomes smaller with increasing temperature and with increasing 
exposition time to an electron or photon beam. Furthermore, it was found that suspended bilayer and 
trilayer graphene present a roughness exponent of £ ~ 0.8. Scanning tunnelling microscopy studies [196] 
further confirmed the static nature of ripples in suspended graphene. The ripples found in Ref. [196] had 
a characteristic lateral size of 5 nm and showed to be static for imaging times of about 5 minutes. The 
used samples were subject to annealing and reported to be free of contamination over areas of tens of nm 2 . 
The amplitude of the corrugation of graphene was further studied in Ref. [177] with transmission electron 
microscopy, where, by analysing the width of diffraction peaks, it was obtained a height and normal variances 

of yj ( h 2 ) = 1.7 A and yj\Vh\ 2 = 0.011, with a characteristic length scale of 10 nm, for suspended samples 
with a size of 0.5-5 pm. It was also found that, at least for length scales smaller than 100 nm, the height 
fluctuation decreases with increasing temperature. 

The exact origin of the static corrugation presented in graphene samples is still unclear. While the 
temperature dependence of the corrugation measured in Refs. [195, 177] indicates that thermal activation 
of phonons should play a role, the theory of flat anharmonic crystalline membranes alone is unable to 
explain the static nature of the observed ripples in [176, 193] and the non-universal roughness exponents 
measured in Ref. [195]. A possible origin of the ripples is that in a suspended graphene sample, its edges are 
constrained and therefore the membrane is unable to freely expand or contract, leading to buckling and the 
formation of ripples [165, 166]. Therefore, residual thermal or external stresses could lead to the formation 
of ripples. Thermal and external stresses have indeed been used to form ripples in graphene suspended over 
trenches as reported in Ref. [10]. This has been explained with the theory from Ref. [197] based on classical 
elasticity, but seems unrelated to the short wavelength corrugations reported in Refs. [176, 193, 195, 177]. 
Another possibility is that in the fabrication process of a suspended graphene sample, relaxation of a 
constrained graphene membrane after removal of the substrate could lead to rippling, or the ripples might 
themselves be inherited from the substrate [198]. Electron-phonon interaction in graphene might also lead 
to a crumpling instability of the membrane [182, 178, 179]. Disorder has also been proposed as a source 
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of ripples. Indeed, the numerical simulations from Ref. [172] show that a 20% adsorption of OH molecules 
leads to the crumpling of a graphene membrane compatible with amplitudes and typical sizes of corrugation 
observed in experiments. Scanning transmission electron microscopy also shows that the presence of point 
defects as vacancies and adatoms, such as hydrogen or carbon atoms, favours rippling [199, 200]. These 
disorder induced ripples might still display dynamics due to migration of the defects. 

Scanning tunnelling microscopy studies were also capable of probing time-dependent height fluctuations 
of suspended graphene [201]. At small bias voltages (~ 0.01 V) and currents (0.2 nA), random on-site 
fluctuations were attributed to thermal fluctuations. These displayed a variance of \J ( h 2 ) = 1.47 nm, with 
a characteristic decay time of r = 8s ({h(t)h( 0)) oc e _t//T ). It was also found that for larger currents and 
gate voltages, effects of thermal expansion due to the local heating of graphene and the electrostatic forces 
induced by the tip can lead to periodic oscillations of the membrane and to mirror buckling [201, 202, 203]. 

3.2.2. Anomalous elasticity 

The Young modulus of graphene was first measured using a blister test in Ref. [25], where an atomic 
force microscope is used to perform a nanoindentation in a graphene membrane suspended over a circular 
hole, with sizes of 1 and 1.5 /im in diameter. The two-dimensional Young modulus, Y e ffi was determined by 
fitting the force-indentation curve by [204, 205, 206, 207] 


F = a 0 n5h + Y eS ^^~, (55) 

where F is the applied force, Sh is the displacement of the membrane at its centre, <jq is a pretension, d is 
the membrane diameter and d = (l.05 — 0.15^ — 0.16z/ 2 ) , where v is the Poisson ratio of the membrane 

[204, 206]. By taking v = 0.165, the obtained value for the Young modulus of graphene was 342 N/m [25], 
a value that makes graphene the stiffest material to have ever been measured. 

Equation (55) is derived from the Foppl-von Karman plate theory, which is described by a free energy 
of the form of Eq. (17) [21]. However, while the high bending rigidity of a plate strongly suppresses thermal 
fluctuations, which can therefore be neglected, on a membrane long wavelength vibrations are thermally 
excited and interact due to strong anharmonic effects. As discussed in Sec. 3.1.1 this leads to an anomalous 
dependence on momentum of the displacement correlation functions, which translate into a scale dependence 
of the elastic constants. As a matter of fact, while the starting potential energy that describes the membrane 
is a local functional of the displacements, the free energy will become non-local [161, 162]. Assuming that the 
strain induced by the indentation is nearly uniform over most of the sample, we can assume that Eq. (55) is 
still valid, but with a scale dependent Young modulus Y e ff(^) oc £~ Vn (Eq. (47)), where l is the characteristic 
length scale of the problem (which can be the size of the membrane, a length scale imposed by the induced 
strain, or a length scale due to disorder or an initial corrugation of the membrane). This justifies the use of 
Eq. (55). 

It has been experimentally found that the Young modulus of graphene increases with the density of 
defects created by bombardment of Ar + ions [208], provided the density is not too high. The measured 
Young modulus, changed from 250 - 360 N/m for pristine graphene, to 550 N/m for graphene with an 
average defect distance of 5 nm (0.2 % coverage of defects). The increase of the Young modulus has been 
explained in terms of the classical theory of flat anharmonic crystalline membranes described in Sec. 3.1.1. 
According to it, the Young modulus becomes smaller at larger length scales. The creation of defects would 
introduce a new length scale, the average distance between defects, which acts as an effective membrane 
size from the point of view of the graphene phonons. Being smaller than the membrane size, this new scale 
would restore the value of graphene’s Young modulus to a higher value than in the absence of defects. The 
Young modulus in Ref. [208] was measured with a nanoindentation technique and by fitting the obtained 
force-displacement curve to Eq. (55) for different graphene drums. In Fig. 3 the experimental value of the 
Young modulus obtained in this way is plotted versus the density of defects for the different drums. The 
Young modulus grows up to twice its pristine value for all the drums with the same law; the experimental 
curves then show a crossover to a decreasing regime as the density of defects is increased, as expected from 
first principle calculations [209, 210, 211]. These competing behaviours can be encoded in a qualitative fitting 
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to the experimental data 


(56) 


Y eS = K (b+ ^ + n^j ^1 - c (J 2 + j 

where rii is the density of defects, K and c are constants to be evaluated from the fit, b is a geometrical 
parameter of the order of the inverse of the area and Iq is the localization length for the flexural phonons 
in pristine graphene. Best fitting with the experimental data, shown in Fig. 3 with the broken green line, is 



Defect density(cm' 2 ) 


Figure 3: Young Modulus as a function of the defects density [208] (see main text). The broken green line represents the best 
fitting to the data via Eq.(56). 

given by the value of the parameters K = 1.5 x 10 9 N m Vu ~ 2 , Iq = 70 nrn, r) u = 0.36 and c = 1.2 x 10 -18 m 2 . 
The localization length Z u takes into account other sources of disorder that can introduce an effective length 
in the problem. The radii used in Ref. [208] vary from 500 nm to 3 p nr and thus the estimated value is such 
that 1 /Iq b. This justifies the independence of the growth of the Young modulus observed. The above 
interpretation is supported by recent experimental measurements of the evolution of the Young modulus with 
external strain [212], where it was found that the in-plane elastic constant increases by a factor of 2 when the 
membranes are under tension. This effect is expected from the elastic theory of membranes, which predict 
that external strain suppresses thermal out-of-plane fluctuations, with the corresponding enhancement of 
the Young modulus [167]. Recent atomistic Monte Carlo simulations have shown that external tension can 
indeed cause an increase of the Young modulus in a factor of 2 [213], in agreement with the experiments [212]. 
However, the discrepancy between the Young modulus measured by indentation experiments on tensioned 
membranes, and that obtained by ab initio calculations is an open question. Another experimental evidence 
of the importance of anharmonicities on the elastic properties of atomically thin membranes has been 
reported by Blees et al. [214], who have shown that the bending rigidity of micro-ribbons of graphene is 
three order of magnitude larger than its microscopic value. These results have been studied theoretically 
in the framework of a renormalization group analysis of flexural phonons on elongated ribbons, which are 
enough to explain the qualitative behaviour [215]. Nevertheless, further studies, both experimental and 
theoretical, are required in order to fully understand the elastic behaviour of graphene. 

3.3. Thermodynamic properties 
3.3.1. Thermal expansion 

Thermal expansion is an intrinsically anharmonic effect. In bulk three-dimensional crystals, thermal 
expansion is well described with the quasi-harmonic approximation. This is essentially a minimal mean 
field extension of the harmonic theory, where the phonon frequencies depend on the lattice expansion. The 
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dependence of the phonon frequencies on the expansion is encoded in the so called Griineisen parameters. 
For crystalline membranes anharmonic effects will be too strong to treat in this matter. At not too high 
temperatures, thermal expansion will be controlled by the low energy out-of-plane phonons and can be 
studied using the elasticity model (17). In a quasi-harmonic approach, the Griineisen parameter for out-of- 
plane phonons is given by [216] 


7(9) 


A + fi 

2 hiq 2 5 


(57) 


which diverges for small momenta. The areal thermal expansion is defined as being the change with temper¬ 
ature of the membrane area projected onto the reference plane. Within elasticity theory the relative change 
of projected area is given by ( diUi ). In terms of the out-of-plane mode Griineisen, the thermal expansion in 
the quasi-harmonic approximation is given by [216, 9] 
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(58) 


From the previous expression, it is clear that in the quasi-harmonic approximation the thermal expansion 
for an infinite membrane will be divergent at any finite temperature [216, 148]. Inclusion of anharmonic 
effects beyond the quasi-harmonic approximation corrects this unphysical result. In the context of the model 
described by Eq. (17), the areal thermal expansion can be exactly written as a normal-normal correlation 
function [216, 27, 148] 

aA = ~ \ Wf ( d M x ) d M x )). (59) 

This is the formal expression of the membrane effect [217]: the intuitive picture that the thermal excitation 
of out-of-plane vibrations should lead to a reduction of the area of the membrane. Equation (59) always 
predicts a negative thermal expansion, being a limitation of the model (17). Higher order anharmonic 
effects, such as cubic terms in diUj, would give origin to competing terms, that at high enough temperature, 
could change the sign of the thermal expansion. Nevertheless, at low temperature, Eq. (59) is the main 
contribution. If the average in Eq. (59) is taken with respect to the harmonic theory, one recovers the result 
from the quasi-harmonic approximation. Going one step beyond the quasi-harmonic approximation, one 
can estimate the contribution due to strong anharmonic effects at low momentum by taking the Ginzburg 
momentum scale (43) as a low momentum cutoff, obtaining [216] 


ola - - —log (A/g*), (60) 
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where A is a large momentum cutoff, which one takes to be the thermal wavelength for flexural phonons 
at intermediate temperatures and the Debye momentum at high temperatures. Using typical values for 
graphene one obtains a thermal expansion of a a ~ —10 -5 K -1 . This estimation is in agreement with 
more involved calculations using a quasi-harmonic approximation based on ab initio calculations [218], 
classical Monte Carlo simulations [26], non-equilibrium Green’s functions method [219], ab initio molecular 
dynamics [12, 220] and self-consistent-field methods [181]. The thermal expansion at very low temperature 
is dominated by flexural phonons, and it was found that in the quantum limit, inclusion of anharmonic 
effects renders the thermal expansion finite in the thermodynamic limit, behaving with temperature as 
a a oc — f or flexural phonons with dispersion relation of the form oc g 2_ W 2 [148, 149]. 

There is a debate regarding the behaviour of the thermal expansion coefficient of graphene with tem¬ 
perature. From the theoretical point of view, while Refs. [218] and [12] predict an always negative thermal 
expansion (Ref. [12] predicts this only for large enough samples, while for small samples there is a change 
of sign), Refs. [26, 219, 220] predict a change of sign in the thermal expansion for temperatures in the range 
600 — 900 K. This possible change of sign occurs due to a competition of thermal contraction due to the 
anharmonic effects involving the out-of-plane mode and thermal expansion due to the in-plane modes. It is 
very important to emphasize that the thermal expansion of graphene greatly depends on whether it is free 
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or if it is supported on a substrate [12, 27, 220]. It was also found that the calculated thermal expansion 
is dependent on the size of the simulated system [12, 220], which can be understood since the anharmonic 
coupling between out-of-plane and in-plane vibrations lead to long range interactions [3] which might be dif¬ 
ficult to simulate numerically. The discrepancy of results can be attributed to the different approximations 
involved in the calculation. An accurate description of the physics would have to take into account quantum 
effects, anharmonic effects in a non-perturbative way and should be extensible to the thermodynamic limit 
(a non-trivial task, from the numerical point of view, since anharmonic effects are stronger at larger length 
scales). 

Experimentally, the thermal expansion coefficient of graphene has been measured in suspended samples 
[10], and was found out to be negative up to 350 K with a value of a.A ~ —7 x 10 -6 K -1 at 300 K. 
The thermal expansion coefficient was also measured in supported samples by Raman spectroscopy [11] and 
was found to be negative in the range 200 — 400 K, with an obtained value of ola ~ —8 x 10 -6 K _1 at 
room temperature. The data from Ref. [10] was successfully described by the theory of Ref. [181], which is 
formulated in the thermodynamic limit, taking into account anharmonic effects in a self-consistent way and 
including lowest order quantum corrections to the classical result. 

The areal thermal expansion discussed so far, being a change of the project area of the membrane, 
is related to the change of the lattice parameter. Another interesting quantity is the change in nearest 
neighbours carbon-carbon distance. There is also no consensus regarding the behaviour of this quantity 
with temperature. Monte Carlo simulations indicate that for small temperatures there is also a reduction 
in the carbon-carbon distance, reaching a minimum at a temperature of T ~ 700 K [26]. On the other 
hand, ab initio molecular dynamics simulations show that the carbon-carbon distance always increases 
with temperature up to T = 1000K, even when the lattice parameter contracts [12]. Experimentally it is 
verified that the carbon-carbon distance increases in graphene supported by an iridium substrate over the 
temperature range of 300 — 1200 K [12]. However, one must point out, that theoretically it is also expected 
that the carbon-carbon distance should increase for supported graphene [12] and therefore, it is not clear if 
the experimental result on supported graphene can be extrapolated to the suspended case. 


3.3.2. Specific heat 

Another thermodynamic quantity that is affected by anharmonic effects is the specific heat. Anharmonic 
effects are significant both in the high and low temperature limits. At high temperature, anharmonic effects 
lead to a violation of the Dulong-Petit law: instead of a saturation of the specific heat to a value of k B per 
degree of freedom, one obtains a linear correction on temperature for high enough temperature. This was 
investigated via Monte Carlo simulations in Ref. [26], where it was found that the specific heat at constant 
volume in graphene per carbon atom for temperatures T > 800 K is approximately given by 

ca = 3 k B ^1 + . (61) 


with Eq = 1.3 eV. Inclusion of quantum effects to lowest order [181] leads to a specific heat of the form 


ca = 3 ks 


k B T 

Ei 



(62) 


with Ei ~ 0.053 eV and E 2 — 0.022 eV. The second term of the previous expression corresponds to the first 
quantum correction to the classical anharmonic result. 

At low temperatures, anharmonic effects also lead to a correction to the power law dependence of the 
specific heat with temperature, since the low energy acoustic phonons acquire a renormalized dispersion 
relation. This was studied in Ref. [148, 149], where it was found that at very low temperatures, where the 
main contribution to the specific heat is given by the flexural phonons, the harmonic result for the specific 
heat, c p oc T, is corrected to c p oc T 4 ^ 4 ~ Vh ), for flexural phonons characterized by a dispersion relation of 
the form ui? oc q 2 ~ Vh / 2 . 
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3-4- Topological defects 

A closely related problem to strain in graphene is the presence of topological defects in the otherwise 
perfect crystalline system. The literature of topological defects in graphene is so vast that one can consider 
that this topic might deserve a review by itself. However, there are several reasons that motivate the 
inclusion of a section about these defects in the present review. If we think on a graphene layer as a 
crystalline membrane, structural defects play an important role in the way a membrane melts. Just as 
elastic deformations affect the electronic properties of graphene at low energies, so will topological defects, 
but in a different way. 

The role of topological defects in the melting and crumpling thermal phase transitions in two dimensional 
membranes was envisaged and developed long time before the appearance of graphene [221, 222]. In the 
melting scenario of Refs. [221, 222] for strictly two dimensional crystals, the positional quasi-long range 
order is destroyed at high enough temperature by the sequential release, first of dislocations and, at higher 
temperature, of disclinations. The possibility of buckling in the out-of-plane direction leads to a reduction 
of the formation energy of defects such as dislocations and disclinations [3, 2, 223]. 

In graphene, due to the covalent nature of the a bonds, responsible of the structural properties, the 
melting transition shows some differences both at qualitative and quantitative levels. Although in graphene 
the presence of topological defects made of a single heptagon-pentagon pair (dislocations) is allowed, it is 
the presence of Stone-Wales defects (formed by two pentagons and two heptagons) that plays a crucial role 
according to Monte Carlo simulations, since they are the structural defects with smallest formation energy 
(in practical terms, the Stone-Wales defect consists in just rotating one carbon-carbon link by 90 degrees). 
The melting transition temperature has been estimated to lie in the range T c ~ 4500 — 5800/Q224, 225]. 
At the transition temperature, the system starts to form an amorphous network of one dimensional carbon 
chains. The same effect has been studied in carbon systems with small number of carbon atoms[226]. We 
should have in mind that the mentioned Monte Carlo simulations are performed keeping constant the number 
of carbon atoms during the simulation. It is quite likely that at lower temperatures than the ones mentioned 
above something more drastic, such as burning, will take place. 

At low temperatures (well below any melting transition) topological defects still have an effect in the 
elastic properties of graphene. Early signatures of the role played by defects on the elastic and plastic 
properties of carbon-based structures appeared in the literature of nanotubes[227, 228]. It was explained 
that the way a carbon nanotube releases tension is to form Stone-Wales defects when the strain exceeds the 
critical value of five per cent. Quite remarkably, the Stone-Wales defects act as nucleation centres for the 
formation of dislocations, in the sense that when more tension is applied, the Stone-Wales defect splits and 
the two pentagon-heptagon pairs start to separate, and even to form more complex defects. 

Back to graphene, it has recently been observed through high-resolution transmission electron microscopy 
that dislocations carry with them sizeable elastic deformations [229]. Interestingly, the dislocations and the 
associated elastic distortion are dynamical, that is, they can move with time in order to release stress. The 
time scales of these movements are of the order of seconds, too slow to have any impact in the electronic 
dynamics, but this defect movement might have an impact in the elastic dynamics of graphene as a mem¬ 
brane. Also, understanding the way these defects evolve constitutes an interesting problem on its own. In 
Ref. [229] it has been observed that, as happens in nanotubes, under stress the plastic deformations start 
with the formation of Stone-Wales defects[230]. Increasing the applied stress, the two pentagon-heptagon 
pairs split apart [231]. In the experiments shown in Ref. [229], the two heptagon-pentagon defects are sep¬ 
arated by several lattice spacings. It has also been found experimentally that the dislocation movement is 
accompanied by the removal of carbon atoms (activation energy of the order of 5 eV) to accommodate the 
energy irradiated by the transmission electron microscope. It has been also observed[232] that buckling has 
a profound impact in the dynamics of a dislocation from its origin, since for a membrane, distortions along 
the third dimension are allowed and graphene samples might also release elastic energy via buckling. This 
buckling process induces a long-ranged interaction between the heptagon-pentagon pairs, and the relative 
direction of the associated Burgers vectors determines the buckling profile, indicating an unexpected large 
feedback between the dislocation configuration, its associated buckling and the way the pentagon-heptagon 
pairs mutually interact through this buckling. 


28 


From the theoretical side, the dislocation dynamics in graphene is still an open problem. Significant 
advances have been achieved with discrete[231] and continuous[233, 234] modelling. Nevertheless, the fact 
that the defect dynamics seems not to be well described by the standard theory of dislocation dynamics 
indicates that the particular lattice topology of the honeycomb lattice might play an important role. 

In the previous discussion, we only gave a small glimpse of the role of small topological defects on 
the elastic properties of a graphene membrane. However, it is by now apparent that in realistic samples 
extended defects are also present and they possibly have a more unexpected impact on the elastic properties 
of graphene. One-dimensional extended topological defects usually appear in chemically vapour deposited 
samples grown over metallic substrates, due to the fact that small graphene grains with different crystalline 
orientation serve as nucleation points and form large-area samples[235, 236, 237]. 

Usually, the structure of such grain boundaries consists on heptagon-pentagon pairs along the grain 
boundary [238]. In Ref. [238] it was stated that if the borders of adjacent graphene flakes are of zig-zag 
type, the Burgers vectors associated to each heptagon-pentagon pair are parallel, while if the borders are of 
armchair type, the associated Burgers vectors lose this constraint. Also, depending on the relative mismatch 
between the crystalline orders, the grain boundaries will be made of well separated heptagon-pentagon 
pairs (low tilt angles) or a dense array of them (high tilt angles). By using molecular dynamics plus ab 
initio calculations, it has theoretically observed that the strength of the graphene sample counterintuitively 
increases with the number of defects in the grain boundary, with the critical strain needed to break the sample 
growing with the number of pentagon-heptagon pairs present in the grain boundaries) [238]. Interestingly, 
calculations based on continuum mechanics fail to explain this behaviour, indicating that it is the microscopic 
lattice structure of these grain boundaries the ultimate responsible of the obtained behaviour. The reason 
behind this is that what controls the failure under strain (which carbon bonds start to break first) is the 
way the strain is not uniformly accommodated along all the carbon bonds, and it tends to concentrate in the 
carbon bonds on the heptagon rings, being then the first to break. This indicates that the system releases 
strain by breaking carbon bonds around the heptagonal defects, leaving the rest of hexagonal rings much 
more free of strain. 

This scenario has been partially questioned due to the fact that graphene can also release strain by 
buckling (or warping) as it has been mentioned previously. Other groups employing similar numerical 
techniques find that actually the stress-strain curve is not monotonic with the increase of defects, meaning 
that not always the strength of the sample increases with the number of defects at the grain boundaries [237]. 
In previous paragraphs we mentioned that the buckling around two close pentagon-heptagon pairs strongly 
depended on the relative orientation of their Burgers vectors, and the same applies here. The relevant 
point here is then to know if the relative orientation between the Burgers vectors plays the same role when 
the distance between heptagon-pentagon pairs is smaller than in the previous cases of a dilute density 
of dislocations. While for zig-zag borders this constraint seems to be quite general, we have much more 
freedom to choose the way the Burgers vectors are not parallel to each other. In this way, for armchair grain 
boundaries, many possible defect configurations with different Burgers vectors can be found, and accepting 
that the buckling profile depends on the Burgers vectors orientation, it is reasonable, at least qualitatively, 
to expect that the way the graphene sample will release strain by buckling (and the behaviour of the strength 
of the graphene sample) is not as uniform and universal as one might initially think. 

4. Electronic properties 

4-1- Strains and the mobility of carriers in graphene 

Graphene is an excellent conductor, with carrier mobilities similar, or higher, than those found in good 
metals, like copper or silver. The lack of dependence of the mobility on carrier concentration in graphene 
strongly limits scattering mechanisms. Defects which induce a short range potential lead to a mobility 
proportional to the inverse of the concentration, which is incompatible with observations. Static strains 
are a possible source of disorder, and limit the mobility of the carriers. Other proposed mechanisms which 
lead to the correct dependence of mobility on carrier concentration are charged impurities in the substrate, 
and resonant scatterers[239]. It is worth noting that thermally excited phonons lead to random strain 
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distributions which can scatter electrons. The in-plane acoustic phonons of graphene are only relevant 
at high temperatures [240], but out-of-plane flexural modes in suspended samples change significantly the 
mobility, being the primary electron scattering mechanism at room temperature [150], giving origin to a 
T 2 dependence of resistivity [241]. In supported samples, out-of-plane phonons become quenched and their 
contribution to resistivity is severely reduced [27]. 


4.1.1. Strains and mobility. 

In order for strains to give the correct (in)dependence of the mobility on carrier concentration, they 
should be correlated over long distances, as initially proposed in [242]. In the following, we analyse the 
influence of strains in electronic transport, following recent experimental and theoretical work reported in 
[243]. As discussed there, it is natural to expect long range correlations in the strain distribution, and 
strains are the only mechanism compatible with weak localization experiments and the observed correlation 
between the mobility and the strength of puddles near the neutrality point. 

The carrier mobility n c is defined as the conductivity per carrier, 

Me = —, (63) 


where a denotes the conductivity, e is the electric charge, and n the carrier density. The highest mobilities 
are found in exfoliated flakes, which are single crystals. They lie in the range 10 3 -10 5 cm 2 /(Vs), and are 
independent of the carrier concentration [244, 245, 246]. Different substrates lead to different mobilities, 
and the highest values are found in samples on boron nitride substrates and in suspended samples, where 
mean free paths can be of the order of few microns. 

If we describe the scattering of carriers by disorder in terms of an effective scattering time r, we can 
write 


e 2ttvft 
^ c h ItF 


(64) 


where h is Planck’s constant, Vp is the Fermi velocity, and kp is the Fermi wavevector. Below, we review 
some properties of the mobility [7, 247, 248], and (following the work in [243]) present arguments which 
suggest that the dominant mechanism which determines the mobility is the presence of random strains in 
the system. 


4-1.2. Form of the disorder potential 

The internal structure of the wavefunctions of carriers in graphene implies that, within states in a given 
valley, an internal degree of freedom, the pseudospin, can be defined. The value of the pseudospin is deter¬ 
mined by the carrier’s momentum, and a variation of the pseudospin affects the phase of the wavefunction in 
a similar way as that of a real spin in a material with strong spin-orbit coupling. The effect of the pseudospin 
can be measured in weak localization experiments, as they probe interference effects by applying a small 
magnetic field which breaks the effective time-reversal symmetry of electronic states around each IF-point. 
These experiments allow for the definition of three different scattering times [249, 250]: (i) To, a scattering 
time which changes momentum but does not influence the phase of the wavefunction, (ii) r*, which describes 
the suppression of interference associated to random fluctuations of the pseudospin, and (iii) Tj, which gives 
the strength of the inter-valley scattering processes (we neglect here the effect of the real spin-orbit coupling 
in graphene, which is very small). The total scattering time, r, which enters in the mobility is given by 


T 1 = T n 1 


T* 1 +T: 1 . 


Weak localization experiments [251, 252, 253] show that r^ 1 ~ r _1 r” 1 . Thus, the random potential 

V ( q) responsible for the scattering processes in graphene must be described by a long range potential, which 
does not allow for inter-valley scattering. In addition, this potential should couple to the pseudospin of the 
wavefunction. 

Further, the experimental observation of a mobility which is independent of carrier concentration implies 
that the effective scattering time t oc kF [see Eq. (64)]. The scattering time is due to processes which 
transfer a carrier from one position on the Fermi surface to another. We assume that the Fermi surface is a 
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circle, and that the scattering mechanism is isotropic. Then, a scattering process by an angle 9 involves a 
momentum transfer q such that |q| = 2kp sin(0), and the inverse of the scattering time (according to Fermi’s 
golden rule) can be written as 

1 _ 2tt N{E f ) 

T H 2 47T 2 

x / dO [1 — cos(0)] (V(q)V(—q))\\^ =2kF sin(e/2) , (65) 

J 0 

where N(E F ) oc k F /v F is the density of states at the Fermi energy E F and () denotes average over disorder. 
The r oc k F dependence mentioned above implies that 

(V(q)V{-q}}\\^ 2 k FSi n i9 /2) oc (66) 

Hence, the potential which scatters the carriers has to diverge like |<f| -1 as k F —» 0. 

4-1-3. Correlation between mobilities and puddles at the Dirac point 

Experiments show a linear correlation between the mobility at high carrier concentration and the min¬ 
imum carrier density which can be measured near the neutrality point [243]. Defects can shift locally the 
position of the Dirac energy, so that, for nominal zero concentration, the carrier density fluctuates, with 
a non zero absolute average. This effect (typically referred to as the existence of puddles [254]) leads to a 
lower bound in the resolution of the carrier density, n*, near neutrality. 

The observed dependence is 1 / p c = Cn*, where C is a constant which does not depend on the type of 
substrate or amount of disorder. 


4.1.4. Carrier scattering by defects and impurities 

Apart from the aforementioned random strains, there are various theoretical models describing the 
charge scattering in graphene. The different microscopic mechanisms relate to different scattering times r 
and carrier mobilities p: 

(i) Weak defects: The simplest scattering mechanism is given by local scatterers which perturb weakly the 
band structure of graphene. Each of these scatterers can be described by a momentum independent 
potential V, resulting in 


t 1 oc N(E F )V 2 oc k F . 


(67) 


(ii) Charged impurities: The electrostatic potential induced by a point charge, 


V(q) = V c (q)/e(q), 


( 68 ) 


where Vc(q) = e 2 /(2eo|<?|) is the 2D Coulomb potential (in SI units), £0 is vacuum permittivity and 
e{q) = 1 + Vc(q)N(E F ) = 1 + (e 2 kp) / (2Trv F So\q\) is the static dielectric function in graphene, satisfies 
the scaling law Eq. (66). The resulting scattering time r and mobility are 


_1 C c hV F 
T OC —;-. 

k F ' 

e 1 


Pc oc 


h c c h 


(69) 


where c c h is the concentration of charged impurities, independent of the carrier density [255, 256, 257]. 
This mechanism can be consistent with the correlation between p c and n* described in Sec. 4.1.3. 
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(iii) Resonant scattering: A strong local potential, like the one induced by a vacancy or an adatom covalently 
bonded to a neighbouring carbon atom, induces a strong resonance near the Dirac energy [258, 259]. A 
potential of this type also satisfies the scaling behaviour in Eq. (66), and leads to a mobility 


H c oc 


e 1 

h Cres 


(70) 


where c res is the concentration of resonant scatterers. 


The scattering time Eq. (67) for weak defects is inconsistent with the observed independence of the mobility 
on carrier density. On the other hand, the long range potential induced by charged impurities, Eq. (68), does 
not modify the pseudospin, while resonant scatterers have dimensions comparable to the lattice spacing, 
and can thus be expected to induce inter-valley scattering. This leaves mechanisms (ii) and (iii) at odds 
with the results of weak localization experiments described in Sec. 4.1.2. 


4-2. Carrier scattering by random strain fluctuations 

The effect of strains on electrons and holes in graphene can be divided into a scalar potential, V s (f) = 
9i {u xx (r) + u vv (r)), and a gauge potential, A el {r) = g 2 [u xx {f) - u yy (r), -2u xy (r)), see Eqs. (19),(20), 
where Uij is the strain tensor, and g\ and g 2 are constants [8] discussed in Sec. 2.3. Random strains will 
give rise to a potential which satisfies the scaling in Eq. (66) if 


{uij(q)Uki{-q)) l k1=fcF oc 


1 

k 2 F 


(71) 


While the gauge field A el does not break the electron-hole symmetry of the graphene bands, so that 
it cannot induce density fluctuations near the neutrality point, the scalar potential V s does break this 
symmetry, and can therefore lead to the formation of puddles [55]. On the other hand, the scalar potential 
does not couple to the pseudospin of the wavefunction, while the gauge potential does, and hence can explain 
the weak localization experiments. The two contributions to the potential lead to the scattering times 
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(73) 


where A[(q) is the component of A el perpendicular to q. as A^ (q) can be gauged away. The average puddle 
density due to the scalar potential is 


1 (Vs{r) 2 ) 1 f ,2 Mg)Vs(-q)) 

7T ( hv F ) 2 47 T 3 h 2 vf J q S 2 {q) 


(74) 


The scalar potential is screened by the static dielectric function, e(q ), while the gauge potential is not. That 
implies that rf 1 <C rf 1 . Hence, when the scattering is induced by strains, r ~ r* ~ r g , in agreement with 
the weak localization experiments described above. 

In the following, we discuss strain due to out-of-plane corrugations and to in-plane displacements. Their 
effect on spectral properties will be discussed in Sec. 5. Both types of strain lead to potentials consistent 
with the scaling behaviour in Eq. (66), and with the observed correlation between mobility and puddle 
density. 
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4-2.1. Random out-of-plane corrugations 

Random out-of-plane corrugations have been observed in graphene on different substrates (where it is 
assumed that they follow the corrugations of the substrate) [245] as well as in suspended samples [176]. The 
fluctuations in the height of the graphene layer, h(f *), give rise to strains Uij oc dihdjh. For the relation (71) 
to hold, we have 

(^dihdjh\f d k hdih\_^ oc t-^. (75) 

This momentum dependence can be obtained if [242] 

(h(q)h(-q)) oc (76) 


The correlation function (76) describes the thermal distribution of classical flexural modes in the harmonic 
approximation (see Section 3.1). Hence, if we assume that the corrugations in a graphene flake are given by 
the quenched shape of that flake at some high temperature, we obtain a mobility which is independent of 
the carrier concentration. 

Assuming a height scaling as in Eq. (76), the relation between mobility and average puddle density 
satisfies 


J_ = n * h_ [ 

H c 71 4e 8e 4 


g|(A + /i) 2 


1 

l°g[l /{k F (n*)a)Y 


(77) 


which, neglecting the last logarithmic factor, only depends on parameters defined in the absence of disorder. 
For reasonable values of g\ and g 2 this relation is in close agreement with the experimentally observed one 
[243], 



Figure 4: (Colour online). Sketch of the potential (left) and forces (right) induced by the substrate on a graphene layer. To 
compare with a floating phase configuration, cf. Fig. 12. 


4-2.2. In-plane random displacements 

The substrate induces in-plane forces on the carbon atoms that constitute the graphene layer, such that 
the atoms will be displaced towards positions which lower the potential energy. This effect is responsible 
for the formation of Moire patterns in graphene on BN substrates [260, 261, 262]. A sketch of a possible 
random potential induced by the substrate, and the resulting displacements, is given in Fig. 4. A detailed 
analysis of this issue concerning periodic deformations is carried out in Sec. 5, where a floating phase is 
assumed. (See in particular the discussion in Section 5.2.) 

Forces F with a modulation K similar to the reciprocal vectors of the graphene lattice, G, lead to smooth 
displacement distributions, modulated by the wavevector k = K — G: 


u{k) oc 


F\\{k) | F ± (k) 
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(78) 
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where Fj|(fc) and F±(k) denote the components of F(k) parallel and normal to fc, respectively, and A and 
p are elastic Lame coefficients, as discussed in Sec. 2.2. The linear strain tensor is given by u^j^(k) = 
i ^ kiUj(k) + kjUi(k)j /2. 

The forces can be written as a function of the potential between the substrate and the graphene layer 

F{K) = iKV R . (79) 


Hence, for K ze G, we find 


U(i,j)(k)u(i tTn )(-k)) oc 
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(80) 


We assume that the potential has random correlations for |fc| <C |G|. Then lim|^T|_^ 0 (Vg_^V_Q + ^j = 
constant, and the correlation between strains follows the scaling in Eq. (71). This mechanism gives a 
mobility which is independent of the carrier concentration. The relation between the mobility and the 
puddles’ density is, similar to that in Eq. (77), 


L = \ ^ 2v f , flf (■, , (A+ 2 fi) 2 
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(81) 


As in Eq. (77), for reasonable values of the parameters of g i and < 72 , the relation between p c and n* agrees 
with experiments [243]. 

Thus, scattering of charge carriers by random strains, as detailed in Sections 4.2.1 and 4.2.2 above, fulfils 
the experimental constraints which a microscopic charge scattering mechanisms in graphene must satisfy. 
Moreover, these constraints, namely (i) the existence of a mobility independent of carrier concentration, 
(ii) the dominance of long range intra-valley scattering processes coupled to the wavefunction’s pseudospin, 
and (iii) the correlation between mobility and puddles’ density near the neutrality point, are not consistent 
with the other models for scattering discussed in Section 4.1.4 [243]. Furthermore, spatially resolved Raman 
measurements show a good correlation between carrier mobilities and amplitudes of strain distributions 
[243]. 


4-3. Electron pumping through mechanical resonance in graphene 

While electron scattering from random microscopic strain configurations or from thermally excited 
phonons contributes to the resistivity of graphene sheets, macroscopic deformations (with amplitudes of 
a few nanometres in samples of micrometres in size) produced in graphene based nano-electromechanical 
systems (GNEMS) can generate electrical current by moving electrons coherently [263]. A typical setting 
for a GNEMS is made up of a graphene flake of length L suspended over a pre-patterned trench. The 
suspended flake serves as a mechanical resonator, which can be actuated via a radiofrequency gate-voltage 
applied between graphene and the substrate [264, 265]. In alternative, optical actuation can be achieved 
by modulating the intensity of a diode laser shone over the flake [264]. Since their first experimental re¬ 
alization [264], GNEMS have been fabricated in various ways: graphene flakes produced by mechanical 
exfoliation[264], CVD[266], or by epitaxial growth on SiC [267] , have been suspended on differently shaped 
substrates. 

In the miniaturization of NEMS to the nanoscale, which calls for low oscillator masses and high resonance 
frequencies [268, 269, 270], graphene with its low mass density, high stiffness, high electronic mobility, and 
its chemical inertness is a promising material [263, 271]. Furthermore, graphene allows tunability of the 
resonance frequency ujq on a wide range by means of strain: Modelling the graphene membrane as a beam, 
e scales with the square root of the applied strain e, 
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where Y denotes the Young modulus of graphene. Eq. (82) takes into account only the effect of tension 
on the fundamental frequency. A realistic model of a graphene resonator must also further consider the 
electrostatic force between graphene and the substrate, the elastic restoring force and possible dissipative 
effects[263, 272], 

Read-out of GNEMS can be performed through electrical or optical means[271]. It has been shown that, 
due to the long wave-length strains induced by vibrations, a GNEM close to resonance can work as an 
adiabatic quantum charge pump [263]. If inversion symmetry is broken by maintaining different chemical 
potentials at graphene contacts on opposite ends of the trench, one obtains a nonzero pumping current such 
that the transport of a single charge per cycle is feasible. Coulomb blockade in such a device would result 
in the transport of an integer number of charges per cycle, so that the ratio between current and frequency 
will be quantized. This paves the way to the application of GNEMS as prototypical standards for current 
and resistance. 

4-4- Other effects of strains on the electronic structure of graphene. 

Strains modify the hoppings between carbon n orbitals. As mentioned in the introduction, this effect 
shifts the Dirac point from the corners of the Brillouin Zone, and leads to the appearance of effective gauge 
fields [8]. When the strains are spatially modulated, an effective magnetic field s induced. A smoothly 
increasing strain leads to an effective constant magnetic field [57]. This effect was confirmed by observing 
the local density of states of highly strained triangular graphene bubbles on platinum[59]. The combination 
of large strains, e ~ 10%, which change over a small scale, £ ~ 15 nm, and a triangular shape favour the 
formation of large effective fields. Numerous experiments have confirmed the existence of significant effective 
magnetic fields in strained graphene, see, for instance[273, 61, 274]. 

4-5. Effects of electron-electron interactions in strained graphene 

The effect of electron-electron interactions in graphene has been a central research area in the field 
since the early days [275]. At half-filling, and unlike in typical metallic systems, Coulomb interactions 
have a long-range character due to the vanishing density of states at the Dirac point, resulting in an 
enhancement of the Fermi velocity at low energies [276, 277, 278]. On the other hand interactions in doped 
graphene are screened, but can nevertheless lead to remarkable phenomena. For instance, a particularly 
appealing scenario occurs when graphene is doped up to the van Hove singularity, where the divergent 
density of states makes it prone to a number of interesting weak coupling instabilities including d-wave 
superconductivity [279, 280, 281, 282, 283, 284]. 

The interplay between strain and interaction effects opens a different avenue to explore and enhance 
novel interaction related phenomena. For instance, by applying uniaxial strain, it is possible to alter the 
electronic spectrum and bring the van Hove singularity to lower energies [89] and rippled disorder can also 
enhance interactions [198]. In what follows, we focus on how electron-electron interactions of different kinds 
can affect strained graphene and the different phenomena they induce. 

4-5.1. Topological phases induced from the interplay between strain and interactions 

Unlike symmetry broken phases, topological states of matter are described and classified by topological 
invariants [106, 107]. This rapidly expanding and multidisciplinary field was boosted by and largely bene¬ 
fited from research in graphene itself [109, 110], cf. Sec. 2.5 . 

An exciting alternative to access topological phases is through the interplay between strain and electron- 
electron interactions. For instance, an early example was motivated by the magnetic catalysis scenario for 
real magnetic fields in graphene [285]. This mechanism, identified first in high-energy physics, refers to the 
enhancement of the dynamical symmetry breaking by an external magnetic held [286]. In the context of 
graphene it is responsible for the opening of a gap catalyzed by the strong magnetic held, which favours 
the electron-hole pairing formation. For pseudo-magnetic helds, it was shown [287, 288, 289] that strained 
graphene is similarly unstable towards an interaction induced gap, but in this case the ground state breaks 
time reversal symmetry and has a hnite Hall conductivity [109]. 

An experimentally plausible scenario emerges when considering strain configurations that result in strong 
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uniform pseudo-magnetic fields [57, 290]. In this set-up, the uniform pseudo-magnetic field splits the spec¬ 
trum into pseudo Landau levels (PLL) 


E n = sgn(n)V2nhv f /1 b 


t 


it' 


^{^\n\~l,rn{z)-,Sg^{n)(j)\ n \ trn {z)) if 71 ± 0, 
(0, ifn = 0. 


(83) 


with l b = y/%/(eB) being the standard magnetic length corresponding to the pseudo-magnetic field strength, 
B , and 4> n ,m are the wave functions of the n-th level of the non-relativistic Landau level problem with angu¬ 
lar momentum eigenvalue to. The states are four-fold (spin and valley) degenerate. They enjoy a Z 2 x SU (2) 
symmetry, where Z 2 represents the time-reversal symmetry while SU(2) are the spin-rotations. Note that 
this is in sharp contrast to the SU (4) symmetry in an external real magnetic field where also the wave 
functions are different for different K and K' points; ip^m 7 ^ i’nm- 


Long-range electronic interactions can lift the large degeneracy of a PLL defined in Eq. (83) and lead to 
interesting ordered states, including topologically non-trivial phases. To tackle this many-body problem, it 
is generally assumed that interactions are not strong enough to mix different PLL. In this case it is possible 
to distinguish the two cases in (83), i.e. n / 0 and n = 0. Considering n/ 0, the projected Hamiltonian 
describing long-range interactions can be written as [291] 

Hi„t = V c(q)ri(q)n{-q). (84) 

Q 

Here Vc(q) is the Fourier transform of the long-range Coulomb interaction and the projection is encoded 
in the definition of the density operator n(q) = rc KjCr (q), where h KjCr (</) are density operators projected 

to the n-th PLL that depend on the valley k = K. K' and spin a =t,4- degrees of freedom. In particular, 
the pseudo magnetic held has opposite signs at the K and K' points, a fact that is reflected through the 
dependence on k. 

Within the Hartree-Fock approximation the leading instability for a 1/4 or 3/4 filled n 7 ^ 0 PLL is a 
state with broken Z 2 valley symmetry. Such a state is not time-reversal symmetric and has a quantized 
&xy = e 2 / h. In turn, for a half-filled n / 0 PLL with spin-orbit coupling a similar analysis shows that 
interactions drive the system to a quantum spin-Hall effect with a quantized spin Hall response [291]. 
Interestingly, the possibility of realizing the elusive fractional Chern insulators (FCI) (for a review see [292]), 
zero-magnetic held analogues of the fractional quantum Hall effect, has also been proposed within this sce¬ 
nario [293, 291]. Via exact diagonalization it was shown that a partially filled n = 0 PLL stabilizes a valley 
polarized fractional quantum Hall state at 2/3 filling through long-range Coulomb interactions that is robust 
to repulsive nearest neighbour interactions [293]. By further tuning the latter to attractive interactions, a 
time-reversal symmetric fractional topological insulator state is stable as well as a spin-triplet supercon¬ 
ducting state. Motivated by these hndings, the stability of the former against different interactions was also 
further addressed [294]. 


Inhomogeneous pseudo-magnetic held configurations can also lead to topological insulating phases. For 
instance, consider a site dependent hopping amplitude of the form [288] 

tij = e x ^te~ x ^ (85) 

with i £ A and j £ B and t a constant hopping amplitude. The function x(77.), with 7Z dehned as a radial 
coordinate from a central hexagon (see Fig. 5) models a lattice gauge potential that can interpolate between 
an approximately uniform pseudo magnetic held (x(77) oc 7 Z 2 ) to a bell-shaped held (x(77) c< In TZ) while 
preserving C 3 rotational symmetry. In analogy to the pristine unstrained graphene [295] fermions at half 
filling subjected to a strain pattern and repulsive interactions can realize topological Chern and spin Hall 
insulator phases. The starting point is the Hamiltonian 

T~L — / ] tij T h.c.^ T U / ] TO,< 7 ^ 2 ,er' T V[ / ] j'(ji “t“ V 2 ^ (^^) 

(i,j) i.er.cr' ({ i,j)),a,a' ((i,j)) 
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Figure 5: Pictorial representations of the hopping amplitudes defined by x(7?.), an increasing of the minimal number of bonds 
1Z to reach a central hexagon. The bond strength is such that if 1Z is odd then x(^) > X(B) e ^ se x(^) = x(-^)j a strain 
configuration that can interpolate between approximately homogeneous to inhomogeneous magnetic fields [288]. Sections 
(a,c,e) and (b,d,f) map into each other under a 2n/3 rotation. 


where tij encodes the strain pattern and { U. V\, Vj } represent the Hubbard on-site, nearest neighbour and 
next-nearest neighbour repulsion strengths respectively. For spinless fermions at half-filling, inducing a 
topological phase can be achieved by developing a complex bond expectation value 

V{(ij» = (4 c j) e C ( 87 ) 

that breaks time-reversal symmetry and acts like an order parameter of the Chern insulator state [109, 295]. 
Being an effective second-neighbour hopping strength, it is the term proportional to V 2 in ( 86 ) that it is ex¬ 
pected to drive the transition. Indeed, under a mean field decomposition of V 2 and with U = V\ = 0, it was 
shown that the strain in Fig. 5 and defined by (87) develops rjij 7 ^ 0 for sufficiently strong V 2 > V c ~ t [288]. 
In the cases with non-uniform strains the order parameter is finite in the vicinity of regions with finite 
pseudo-magnetic field flux. 

Spinful fermions allow for other competing orders to emerge, including the spin Hall effect. In this case the 
on-site Hubbard interaction U favours a spin-polarized ferromagnetic state for small V 2 [293, 296] that can 
turn, upon increasing V 2 , into spin-Hall effect [288, 291]. A nearest neighbour interaction V±, on the other 
hand, favours a charge density wave order which restricts the anomalous Hall phase to appear only when U is 
sufficiently small and V 2 > V\ [296], a situation resembling unstrained half-filled graphene [295, 297, 298]. An 
interesting open question is whether these topological phases survive to quantum fluctuations [299, 300, 301]. 

It is worth noting that the effect of including both the pseudo (b) and external (B) uniform magnetic 
fields on the integer quantum Hall effect in graphene has also been addressed, including the effect of on site 
U and nearest neighbour V] interactions [302, 303]. For B > b and at half-filling, the Freeman splitting 
due to the real magnetic field will lift the spin degeneracy while interactions are left with the task of lifting 
the valley degeneracy (either forming a charge density wave or ferrimagnetic order) triggered by the pseudo 
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flux [287]. The physical effect is a standard quantum Hall conductivity sequence a xy = ne 2 /h, in sharp 
contrast with the non-interacting result a xy = (4 n + 2)e 2 /h with n £ Z that is a consequence of the valley 
and spin symmetry of non-interacting electrons. In the strain dominated regime (b > B) an alternating 
sequence between a xy = 0 and cr xy = =f2 e 2 /h occurs as a function of electron or hole doping. In this case, 
and similar to the B = 0 case discussed above, short range interactions can stabilize quantum spin Hall and 
anomalous Hall insulators among other ordered phases depending on their relative strengths [303]. 


4-5.2. Interaction induced magnetic and non-magnetic phenomena in strained graphene 

Competing with the topological states described in Sec. 4.5.1, other important interaction induced orders 
can emerge in different strained graphene configurations. 

An important case is that of uniaxial strain [90] . Moderate uniaxial strain is known to deform the low energy 
graphene spectrum into an anisotropic Dirac cone resulting in an elliptical Fermi surface [8]. The effect of 
long-range Coulomb interactions in anisotropic Dirac models has been studied (prior to the experimental 
realization of graphene) in QED 3 theories relevant for cuprate superconductors [304, 305] , and was revisited 
recently within the strain context [306]. The Dirac anisotropy can be parameterized by the ratio of the 
Fermi velocities along the x and y directions 


v y /v x = 1 + 5. 


( 88 ) 


Within a renormalization group approach, the anisotropy is found to be irrelevant in the renormalization 
group sense if <5 <C 1, i.e. the isotropic Dirac cone is a stable fixed point [304, 305, 306]. When <5 is sufficiently 
strong a quasi one-dimensional excitonic insulator with 

(tp\ipA - V’sV’s) 7 ^ 0 (89) 


appears to be the leading instability [306]. Furthermore, uniaxial strain can also affect the onset of the 
ferromagnetic instability of neutral unstrained graphene under the long range Coulomb interaction [307], 
driving the instability to smaller couplings within the range of those expected for suspended graphene [308] . 

Due to its possible experimental relevance, it is important to note that uniaxial strain also alters the 
screening properties of graphene that are encoded in the dielectric function e(q, ho) at finite frequency (w) 
and momentum ( q ). In the random phase approximation it is given by 

£{q,iui) = 1 - V c (q)B w (q, iu) (90) 

Here Vc{ q) oc 1/q is the bare Coulomb interaction and nW(q, ilo) is the one loop polarization function. For 
all (q, iio), the latter quantity can be read off from the result in massless quantum electrodynamics (c.f. Eq. 
(2.14) in [275]). The effect of uniaxial strain is apparent by restoring the dependence on the anisotropic 
Fermi velocities 




N v 2 x q 2 x + v 2 v q 2 y 
v yjv x q x +v%q y + ui 2 


(91) 


where TV = 4 is the number of fermion flavours (one for each spin and valley). The effect of such strain 
on e(q, ico) permeates into the behaviour of the van der Waals forces between strained graphene sheets, 
that can be shown to increase with increasing strain anisotropy [309]. Through a renormalization group 
calculation it was shown that, for large separations between undoped strained graphene sheets, the inclusion 
of electron-electron interactions reduces the van der Waals interaction between graphene sheets [309]. 


For interacting graphene in a uniform pseudo-magnetic field, the different ordered phases resulting from 
Coulomb long-range interaction, on-site U, nearest neighbour V\ and next nearest neighbour interaction Vj 
have been extensively studied [293, 296, 288, 291, 310] and discussed in Sec. 4.5.1. From first principles 
however, [311] it is estimated that U is the largest energy scale, which, if strong enough would in principle 
favour a spin-ferromagnet [293, 291, 296]. However, with mean-field theory and Monte-Carlo simulations 
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it was found that, for finite samples, the system develops a global anti-ferromagnet with zero magnetiza¬ 
tion [310]. Although locally this state is ferromagnetically ordered, it displays a magnetization that changes 
sign between the bulk and the edges precisely to deliver a zero total magnetization. This is ultimately tied 
to the fact that although in an infinite system the support of the zero-modes is restricted to one sub-lattice 
(say A), in a finite system there is support for these also in sub-lattice B but only at the edges. The mag¬ 
netization, being tied to the zero modes, changes sign as the support of the modes changes from A (bulk) 
to B (edge) keeping the total magnetization zero. 

To conclude this part we note that long-range Coulomb interactions can reduce the stability of edge states 
resulting from strain [312]. On the other hand an on-site repulsion together with uniaxial strain can lead to 
edge-ferromagnetism [313]. 

4-5.3. Superconductivity in interacting strained graphene 

We finish this section devoted to interactions with a note regarding superconducting states in strained 
graphene induced by attractive interactions. As mentioned above, a spin-triplet superconducting state was 
shown to exist via exact diagonalization in strained graphene [293] with a finite superfluid density 

_ 1 d 2 E g 
ns 2 86 2 ’ 

calculated from the variation of the energy of the ground-state E g upon twisting the boundary conditions 
an angle 9 [314], Moreover, strain can lead to interesting superconducting states in graphene [315, 316]. 
In particular, an exotic / + is superconducting state that spontaneously breaks time-reversal symmetry is 
expected to occur when attractive on site U and second nearest neighbours interactions Vi are comparable 
(JJ ~ Vi) [316]. This state ultimately originates from the competition between the U favoured singlet s-wave 
pairing and the Vi favoured triplet /-wave pairing. 

4-6. Optical properties of strained graphene 

The analysis of light scattered inelastically by graphene by means of Raman spectroscopy has long since 
become one of the standard characterization techniques in graphene research [317]. We give a brief account 
of the effect of strain on graphene’s Raman spectrum in Sec. 4.6.1. However - and not surprisingly for 
an atomically thin structure - the interaction between graphene and light is weak within a broad band of 
frequencies, and elastic scattering is the by far dominant process underlying most of the potential applica¬ 
tions of graphene in optical devices [14, 16, 318, 319]. Thus the optical properties of graphene are mainly 
determined by its electronic structure, which is sensitive to strain. In particular, it has been proposed to 
use uniaxial strain to break the optical anisotropy of two-dimensional graphene, thus converting it into a 
dichroic material (see Sec. 4.6.2). Finally, graphene’s magneto-optical properties, of which, in particular, a 
large Faraday effect has attracted attention recently [15, 16], are determined by the nonlinear Landau level 
spectrum of massless Dirac Fermions. As a zero-field quantum Hall effect and Landau level formation in 
graphene can be realized with strain-induced pseudomagnetic fields [57, 59], the interplay between real and 
pseudomagnetic fields in graphene based optical devices seems to offer interesting possibilities for future 
applications, on which we comment in Sec. 4.6.3. 

4-6.1. Strain monitoring with Raman spectroscopy 

Raman spectroscopy typically uses light in the infrared (IR) or near ultraviolet (UV) spectral range. 
Because the laser energy is large compared to the phonon energy, the scattering mechanism involves electronic 
excitations in intermediate states, rather than direct photon-phonon coupling. Thus, apart from providing 
information about the phonon spectrum, Raman spectroscopy can also shed light on the behaviour of 
electrons [320] and complement transport measurements. As tensile strain usually results in phonon softening 
(and the opposite for compressive strain), Raman spectroscopy provides a tool for strain monitoring. 

The Raman spectrum of defect-free single layer graphene is characterized by two peaks: The G peak, 
which corresponds to the excitation of a phonon in the doubly degenerate in-plane optical mode Ei (see 
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Table 2), and the 2D peak, due to a two-phonon process which involves scattering of electron-hole pairs 
between neighbouring Dirac cones [321]. 

Under uniaxial strain, graphene’s E 2 mode splits into two components Ef 1 , which are perpendicular and 
parallel to the strain axis, respectively. As the E£ mode experiences less softening due to tensile strain, 
the Raman G peak is seen to split into a G + and a G _ band [322, 317]. Both redshift with increasing 
strain, and their splitting increases. Polarized measurements of the G ± intensities allow the determination 
of the crystallographic orientation with respect to a known strain axis. If the uniaxial strain is applied along 
zigzag or armchair directions, the 2D peak likewise splits into two distinct submodes [320]. Magnitude and 
polarization dependence of the mode splitting was found to be consistent with the strain-induced motion of 
the Dirac cones predicted by tight binding calculations [90]. 

In surface enhanced Raman spectroscopy, photonic structures supporting plasmon resonances (such as 
arrays of metallic nanoparticles), are deposited on graphene. The strong electromagnetic nearfield associated 
with the plasmon can increase the intensities of graphene’s Raman spectrum by several orders of magnitude 

[323] . In turn, the Raman spectrum of graphene serves as a detection channel for plasmonic field enhance¬ 
ment. With graphene placed on top of a photonic cavity structure, the slight sagging in the suspended areas 
locally strains the graphene membrane. The resulting phonon mode softening allows to distinguish between 
enhanced Raman signals originating from strained and unstrained areas. Thus the strained graphene mem¬ 
brane can be used as a probe that locally resolves the plasmonic field enhancement of photonic structures 

[324] , 


4-6.2. Optical properties of uniaxially strained graphene 

Elastic (Rayleigh) scattering of light by graphene can be described by the Fresnel equations for a thin 
film with an appropriate sheet conductivity ct(uj) [325, 326]. For normal incident light of frequency u>, this 
yields 


A( W ) = l-|t( W )| 2 -|r( W )| 2 


4cr(o;)/(g 0 c) 

[2 + cr(w)/(g 0 c)] 2 


(92) 


for the absorbance A(u>) of free standing graphene, where r and t denote the Fresnel coefficients for reflection 
and transmission, respectively, c is the speed of light and £0 the vacuum permittivity. 

Modelling the graphene charge carriers as non-interacting, massless fermions in the Dirac-cone approx¬ 
imation yields a universal conductivity of a = e 2 /(4 K) for all ui > 2 Ep/h [327, 328]. When inserted into 
Eq. (92), this results in A = na + 0(a 2 ) ~ 2.3 (where a = e 2 /(Ans 0 hc) is the fine-structure constant) which 
corresponds to the experimental value observed at wavelengths between 450 ran and 4.2 /jm [329, 330]. At 
energies below 2Ep, graphene supports collective plasmon oscillations, which can be described by treating 
electron-electron interactions within the random phase approximation (RPA) [331]. In the long wavelength 
limit hui < 2 Ep, the plasmon dispersion reads 


w(jfc) = (2 ackE F /h ) 1/2 


(93) 


[332, 333, 13]. For typical graphene-on-substrate samples, Ep is of the order of 0.5 eV or lower. This 
fixes the application range of graphene plasmonics to the THz and mid-IR band (wavelengths > 1/mr), 
while graphene based optical broadband applications can operate in in the near-IR and most of the visible 
spectrum (wavelengths between 450 and 1 pm) [14]. 

The effect of uniaxial strain on graphene’s optical conductivity can be described by introducing strain 
dependent hopping parameters into the standard tight binding Hamiltonian [91, 90], which deforms the 
Fermi surface into an ellipse (see also Sec. 2.3). This defines a fast and a slow optical axis, with the latter 
(and, correspondingly, a lowered conductivity) oriented closely along the direction of strain. The frequency 
independence of absorption in the near-IR and visible range stays intact under uniaxial strain, but the 
absorbance of linearly polarized light now depends on the polarization direction 4 >j with respect to the slow 
optical axis [89], 


A ss na[A — 2a\5kp>\ cos 2 <f>j\. 


(94) 
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and the degree of polarization dependence is determined by the strain-induced shift 5kr> = eA s (l + v)/2 
of the Dirac points from their position at K , K' in the Brillouin zone. Here e denotes the magnitude of 
applied strain, v w 0.16 is graphene’s Poisson ratio, and aX s ss 3-4 [91]. The periodic modulation (94), 
measured in the transmission of visible light through bended graphene [334] , allows for an optical detection 
of magnitude and orientation of uniaxial strain, with possible applications in passive, atomically thin, and 
flexible strain sensors. A reduced optical absorption along the direction of strain was likewise observed in the 
Drude response of uniaxially strained graphene at far-IR wavelengths of A = 30-300 /im [335] . Importantly, 
the conductivity modulations resulting from tensile strain up to nearly 20% have been shown to be reversible 
[336, 334], 

The square root character of the plasmon branch [Eq. (93)] as well as its n 1 / 4 dependence on the charge 
density are not changed by uniaxial strain. However, corresponding to the strain dependence obtained for 
the optical conductivity and the resulting absorption in Eq. (94), the plasmon dispersion gets steeper for 
k perpendicular to the direction of strain and is flattened for wavevectors along the strain direction, with 
the amplitude of the modification proportional to e [337]. In a similar fashion, the effect of scalar and 
vector potentials with one-dimensional periodicity [that is, V s = V s (x), A el = A et (x)], induced by periodic 
strain waves, was found to render the low-energy part of the plasmon dispersion anisotropic [338]. Further, 
uniaxial strain (via its effect on the electron-electron interaction) is expected to modify the electron-plasmon 
scattering structures which are visible near the Dirac points in ARPES studies [339]. 

For realistic values of e < 20%, and in contrast to non-uniform strain configurations (as discussed in 
Secs. 2.3 and 5), uniaxial strain is not expected to open a gap (resulting in a vanishing optical conductivity) 
in extended two-dimensional graphene [90]. However, tight binding calculations for graphene nanoribbons 
predict that realistically small uniaxial strain results in a periodic modulation of the bandgap in armchair 
ribbons, while in zigzag ribbons it creates a band gap proportional to the magnitude of applied tensile strain 
[340]. As the band-structure is modified, the optical transition energies between maxima of occupied bands 
and minima of unoccupied bands change, resulting in a shift of the maxima in the optical loss-function 
Z(co) = [Ime(u)] -1 [341]. 

4-6.3. Magneto-optics with strained graphene 

As in the magnetically unbiased case, the transmission of light through a graphene layer in the presence 
of a uniform external magnetic field can be described within the Fresnel formalism [342], with the Fresnel 
coefficients for reflected and transmitted amplitudes now being functions of both graphene’s diagonal and 
Hall conductivity. Analytical expressions for a xy and a xx as functions of frequency and magnetic field have 
been derived in Refs. [343, 344]. Due to the action of the magnetic field on the conduction electrons, an 
incident beam of linear polarization along the x-direction will, after passing the graphene layer, consist of 
both x and y components. The total transmission T through biased graphene is a sum of the intensities 
in both polarizations, and for an incident beam oriented normal to the surface of suspended graphene, we 
arrive at [345, 16] 


T(u) = | t xx \ 2 + \t xy \ 2 = 1 - Regxx{u}) + 0(a 2 ). (95) 

£qC 

Due to the non-equidistant spacing of graphene’s Landau levels, the transmission spectrum Eq. (95), 
as observed by IR spectroscopy [346, 347], displays resonant dips at frequencies where the incident light 
can excite electronic transitions between Landau levels E n = ±vf V2 eHBn (where n = 0, 1 , 2, ...). Calcu¬ 
lations using density-functional theory based tight-binding theory (DFTB) predict that the same Landau 
level spectrum can be obtained at zero external field by twisting graphene nanoribbons with several hun¬ 
dreds of nanometres width [348]. The strain produced by twisting the ribbon around its axis results in a 
pseudomagnetic field B= —Bf ~ y 2 x (where 7 denotes the twist rate, and |x| < W/2 the coordinate 
perpendicular to the ribbon axis) that changes sign when going from one ribbon edge to the other. The 
resulting electronic DOS for zig-zag as well as armchair ribbons is predicted to mimic Landau quantization 
of Dirac fermions in magnetic fields up to 160 T. 
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For graphene, the Faraday effect (and the related Kerr effect), in which linearly polarized light trans¬ 
mitted through (or reflected by) a material in the presence of a magnetic field has its polarization plane 
rotated, is surprisingly strong: monolayer graphene in modest fields has been found to turn the polarization 
by several degrees [15, 16, 345, 342]. With the incident beam as well as the external magnetic field oriented 
perpendicular to the surface of suspended graphene, we arrive at a Faraday rotation angle [15, 345] 

e F (u) = _ 1 Reg *» +0 {a 2 ). (96) 

2 £ 0 C 

dp is large for frequencies where u xy shows resonances. In the classical cyclotron resonance regime, this 
occurs for to ~ Vp\eB / Ep\, while in the quantum Hall regime, laser frequencies matching the transition 
energies between two Landau levels lead to large Faraday angles [15]. 

The Faraday and Kerr effects are of technological importance in the development of optical diodes 
and other non-reciprocal optical elements [349]. For these magneto-optical applications, graphene offers 
a great potential because Landau-level formation and a nondiagonal optical conductivity in graphene can 
be obtained either by applying an external magnetic field or through strain-induced pseudomagnetic fields 
[350, 303]. With an additional strain-induced pseudomagnetic field B s: electrons in the I\ and K' valleys 
are exposed to an effective field B e g = B ± B Sl respectively [350, 303]. An effect of pseudomagnetic fields 
on the Faraday rotation is thus possible at frequencies where the contributions to a xy resulting from the 
two inequivalent valleys do not cancel out in Eq. (96). For hu> <C Ep , the overall dependence of a xy (uj) on 
the effective held scales as B~g [343, 344], which opens possibilities for a strain-induced Faraday rotation in 
the THz band. 


5. Superlattices 

The properties of misaligned graphene multilayers and graphene deposited on other two-dimensional 
crystals have been the subject of intense research recently. The simplest of these, twisted bilayer graphene, 
exhibits different stackings (local atomic alignments between layers) when the two layers are rotated relative 
to each other by a finite angle 9 [351, 352, 353]. Such a rotation gives rise to a periodic spatial modulation 
of the stacking, known as a Moire pattern or stacking superlattice, with a period that goes as ~ a/9 at 
small 9. All sorts of interesting electronic reconstructions arise from the stacking modulation, including low- 
energy van-Hove singularities, secondary Dirac points away from charge neutrality and Hofstadter butterflies 
[354, 355, 356, 357, 358, 359]. These reconstructions can be understood as the consequence of non-Abelian 
gauge fields that emerge in the low energy sector of the system [360], similarly to the Abelian pseudogauge 
fields associated to strains in a graphene monolayer. 

Together with the modulated stacking, a rather strong modulation of the local inter-layer adhesion energy 
density develops in the twisted bilayer. The inhomogeneous adhesion creates a surprisingly strong strain 
modulation in the system, whereby the areas with the preferred stackings (AB, or Bernal stacking, with 
higher adhesion) become expanded at the expense of other regions. This transforms the Moire patterns into 
a characteristic hexagonal mesh of stacking solitons (linear stacking faults, around ~ 10 nm wide), usually 
with the same period as the original (unrelaxed) Moire [351]. Once more, a range of intriguing electronic 
effects emerge from the soliton strain fields, e.g. the development of a helical network of states carrying 
valley currents along the solitons when the bilayer is placed under a perpendicular electric field [361], or the 
suppression of conductance close to neutrality due to scattering on the solitons [362, 363]. 

More complicated Moire structures and inhomogeneous adhesion-induced strains are possible, and even 
ubiquitous in more complex multilayers [355, 208]. Graphene trilayers, for example, can exhibit ABA and 
ABC stacking domains, visible with Raman spectroscopy [364], and separated by compressed stacking soli¬ 
tons. Unlike in the bilayer, these two types of stackings, despite being very close in energy, are electronically 
very different. In particular, ABA domains remain metallic under large out-of-plane electric fields, while 
ABC domains develop a sizeable gap. This enables all-electric manipulation of stacking domains in the 
trilayer, and of their elastic soliton boundaries, using a biased STM tip [365]. These examples clearly show 
that the spatial modulation of stacking alignment in graphene multilayers creates a fascinating playground 
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where elasticity and electronic structure become strongly coupled, and where quite a few problems remain 
open [366, 367]. 

But stacking Moire superlattices are not restricted to graphene multilayers. A simple graphene monolayer 
grown or deposited over different insulating crystals such as SiC and hBN also exhibits Moire patterns, 
stacking solitons and a remarkable electronic structure. High quality samples with outstanding transport 
properties are provided by these setups (cf. Sec. 4.1), and other heterostructures, with [17, 368] or without 
[369] graphene as one of their building blocks. Due to their current interest, we will now review in more 
detail the physics of graphene/hBN structures, as a representative of this class of hybrid Moire superlattices. 

Superlattices arise in this context due to the superposition of 2D crystals with different periodicities, 
stackings and/or alignments. To study their electronic band structure, at least two aspects must be con¬ 
sidered. First, the deformation of the layers’ pristine periodicity that may take place spontaneously due to 
Van der Waals adhesion potentials, which drives the multilayer into a non-trivial equilibrium configuration. 
Second, the electronic coupling between layers in this new scenario, usually regarding the derivation of an 
effective description. 

The Frenkel-Kontorova model addresses the first issue, namely to study the ground state of a system 
where a substrate potential and the elastic energy of a ID or 2D lattice compete. [370, 371, 372, 22] It exhibits 
three possible ground state phases, depending on the relative periodicities of the lattices and the ratio of 
adhesion versus elastic energy: (i) a generalized floating phase, where there are only local deformations due 
to a weak interaction with the substrate, the Moire pattern preserving its undeformed periodicity; (ii) a 
commensurate phase, where there is a modified overall periodicity of the Moire pattern, resulting from a 
finite global deformation; (iii) an incommensurate phase, characterized by the presence of aperiodic solitons, 
and separating two commensurate phases as a function of lattice mismatch. 

These three phases have been experimentally encountered in a variety of heterostructures containing 
graphene. For a twisted bilayer, Ref. [351] observed, in different regions of the sample, both commensurate 
structures with hexagonal symmetry and the structural solitons. For graphene over hBN, Ref. [260] reported 
a phase transition between a commensurate structure and a floating phase as the Moire period exceeds 
~ 10 mn. 3 Fig. 6 is a reproduction of this experimental evidence, with the Young modulus FWHM over the 
Moire wavelength as the order parameter. 

Determining the strain field that carries the system to its ground state is the first issue to be addressed. 
The application of the Frenkel-Kontorova model to two dimensions is not trivial, since it demands the use of 
a complex mathematical machinery, [371, 370] as well as an accurate knowledge of the shape of the adhesion 
potential. To predict theoretically the deformation in this framework is, therefore, a difficult task. 

A first strategy to tackle this problem consists on focusing solely on floating-phase deformations, namely 
strain superlattices with the same periodicity as the Moire pattern. We consider here the simple case of a 
bilayer. Then, such deformations conform minimal commensurate configurations with reciprocal superlattice 
generating vectors[373, 374, 375, 376] 

Gj = 9j - g'j- (97) 

Vectors gj and g'j are the reciprocal lattice generating vectors of each of the two crystals, chosen such that 
the \Gj | are minimal. A sketch of this construction for graphene on hBN (honeycomb lattices with a finite 
mismatch and a relative rotation) is shown in Fig. 7. Gj are often referred to as “first star” vectors when 
dealing with hexagonal symmetry. 

The reciprocal vectors Gj and the corresponding superlattice vectors A :/ define a superlattice that may 
or may not be commensurate with each of the two layers. In the incommensurate case, no true periodicity 
of the bilayer exists, and the basic tenet of Bloch’s theorem cannot be invoked when considering electronic 
structure. It has been shown, however, that for a small value of the graphene lattice constant a over the 
interlayer distance d , incommensurability and exact commensurability yield the same observable electronic 


! Tlie notation in Ref. [260] differs from the customary denominations in the Frenkel-Kontorova framework: they refer to 
the floating phase as incommensurate, in the sense that the top layer does not conform to the substrate. They disregard what 
in this review we call incommensurate phases, since no aperiodic solitons are observed. 


43 





Figure 6: (a) and (b): the Young modulus distribution, measured by AFM, for structures with 8 and 14nm Moire patterns, 
respectively. Scale bars are 10 nm long, (c): Cross-sections of the Young modulus distribution taken along the dashed lines 
in (a) (top) and (b) (bottom), (d): Ratio between the full width half maximum (FWHM) of the peak in the Young modulus 
distribution (as marked by arrows in (c)) and the period of the Moire structure L, as a function of the period of the Moire 
structure for several samples. Reproduced from Ref. [260] with permission of the authors. 



Figure 7: Shortest reciprocal lattice generating vectors of two misaligned honeycomb lattices (red and blue) and of the resulting 
superlattice (green). 

properties. [377] It is then legitimate to approximate any system with a/d —> 0 as a commensurate structure 
with a reciprocal superlattice generated by Eq. (97). This is known as the continuum approximation, 
and allows one to consider the concept of a low energy band structure assuming an effectively periodic 
superlattice, even though Bloch’s theorem does not apply at the smallest length scales. Considering in-plane 
deformations u(f) with the same periodicity (floating phase approximation, with purely local deformations), 
one arrives at the simplest nontrivial description of the coupled electronic and elastic properties of these 
Moire crystals. In subsequent sections we focus on the continuum approximation approach, which turns out 
to be enough to describe a number of systems, such as graphene on hBN. 
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Figure 8: Renormalization of the Fermi velocity in bilayer graphene for small twist angles 8 (strong coupling). a± is given in 
the main text. The inset depicts vf/vF in the weak coupling regime. Red is respective to Eq. (99), whereas the blue curve 
corresponds to more accurate numerical calculations. Reproduced from Ref. [373] with permission of the authors. 


5.1. Undeformed Moire superlattices 

The first studies of graphene superlattices were restricted to floating phases with zero deformations, very 
particularly for the case of twisted graphene bilayers. An effective continuum model for unstrained twisted 
bilayer graphene can be found in Refs. [373, 378, 379]. The Dirac cone structure is preserved, although 
with a renormalized Fermi velocity v* F . For large enough twist angles 0, and at low enough energies, a Dirac 
model remains valid for carriers in each of the two monolayers, which behave as if decoupled, 

'H 0 s = v* F '^2'iljia-kijj^, (98) 

k 


Here ipt. = (at, it), while at ( bt ) are the creation operator of the k Bloch state in the A (B) sublattice. The 
correction to the Fermi velocity v* f /vf as a function of the twist angle is plotted in Fig. 8. In the weak 
coupling regime (large relative rotation angles 9 > 1°), 


V* F ~ VF 


1 — 3 a[y 
1 + 6 a ’ 


(99) 


where a± = kfj , t± is the nearest neighbour interlayer hopping (t± ~ 330meV), kg = 2A"sin(0/2) and 
K = 47r/3a is the modulus of the AT-point wavevector.[373] At energies around v* F kg/ 2, the decoupled 
monolayer description breaks down. Carriers from different layers hybridize, giving rise to a 0-dependent 
van-Hove singularity in the density of states. In contrast with unrotated bilayer graphene, a gap does not 
open in its twisted counterpart when applying a potential difference between the layers. [379] 

As the relative rotation angle 9 falls below 1°, the system enters a more complicated strong coupling 
regime. When neglecting deformations, Dirac cones survive in this regime. However, the corresponding 
renormalized Fermi velocity becomes much smaller and is even completely suppressed for a set of ‘magic’ 
rotation angles, implying a strong enhancement of the density of states at the neutrality point. Ref. [378] 
also points out the possibility of annihilation and regeneration of Dirac nodes based on the ratios of the 
Slonczewski-Weiss-McClure parameters, as well as encountering a quadratic dispersion relation near the 
Fermi energy. 

As for graphene over hBN, a floating phase analysis was conducted in Refs. [380, 381]. Not including in 
plane deformations, ab initio calculations suggest a modulation in the distance between layer and substrate 
(Fig. 9), which in turn results in a gap opening. Interaction between electrons has been shown to enhance 
this gap to the extent of experimentally measured values (~ 30 meV).[380, 382, 383] If, on the other hand, 
the interlayer distance is assumed constant throughout the system and interactions are neglected, the gap 
averages to zero without deformations. 
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Figure 9: Modulation of the graphene distance to hBN as a function of the stacking according to ab initio calculations. With 
AB (BA) stacking, we refer to one sublattice of graphene overlapping exactly with boron (nitrogen) atoms. Reproduced from 
Ref. [380] with permission of the authors. 


5.2. Spontaneous deformations in the continuum approximation 

Strains, such as the spontaneous deformations expected from adhesion interactions with a crystalline 
substrate, induce significant changes in the electronic structure of graphene. [59, 57] Let us assume a certain 
purely in-plane deformation, described at f by an in-plane displacement u(r) = (u x (f), u y (r)) with respect 
to its unstrained configuration. The change in hopping integrals arising from said deformations can be recast 
into the scalar (V s ) and pseudogauge ( A el ) potentials [36] discussed in Sec. 2.3, with the symmetry described 
by Eq. (8) and yielding the Hamiltonian of Eq. (25). 

In the continuum approximation, the influence of the substrate on the electronic structure acts via two 
mechanisms. First, as pseudopotentials ( V s ,A el ) associated to strains, induced in turn by the adhesion 
potential. Second, through the electronic coupling with the substrate. To obtain an effective model out of 
these ingredients one must add the pseudopotentials of Eqs. (19), (20) to a pristine layer, and also integrate 
out the electronic degrees of freedom of its substrate, which creates a self-energy in the form of a local 
potential in graphene for strongly insulating substrates. 

We will focus on the case of small rotation angles between graphene and a hBN substrate, as it has 
caught the attention of most part of the recent literature due to the high quality of the ensemble. Since 
in the continuum approximation we may consider the superlattice to be minimally commensurate for any 
rotation angle, the deformation field u[f) will be assumed to have the same periodicity as the graphene/hBN 
Moire, encoded in Gj. Furthermore, the bulk substrate will be treated henceforth as undeformable (hBN is 
assumed to be a thick and rigid crystal). The smoothest possible deformation field u(r) in the superlattice 
unitary cell thus reads 


u{r) = Y J Z S . 

3 


e iG r r. 


( 100 ) 


Since u(r) is real, its harmonics satisfy ug = u* - . Moreover, the symmetry point group S of the super¬ 
lattice imposes some more constraints, namely Sug = Uq-iq . V«S £ S. For example, when dealing with 

S = C 3 , only the vector (Gi belonging to the superlattice basis) will be independent in Eq. (100). 

Symmetry-based analysis .— Under the foregoing considerations, effective theories dependent on the de¬ 
formation u(r) have been formulated. A first foray into the computation of the Moire band structure in the 
continuum approximation, including deformations, was based on symmetry principles. [375] 

For monolayer graphene over a hBN substrate, the Dirac cone at the Fermi energy is preserved. Moreover, 
due to the folding of the single layer bands over the superlattice 1st Brillouin zone, minibands were predicted 
to exhibit a combination of three possible features: (i) three mini Dirac points (rnDPs) in the conduction 
or valence band, located at the superlattice 1st Brillouin zone edge with momentum HufCj/ 2; (ii) a single 
mDP at the ±K super point (which is analogous to graphene’s K point but respective to the superlattice); 
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Figure 10: Band structure for different choices of the parameters calculated in Ref. [375]. Red, blue and magenta arrows point 
to multiple, single mDPs and triple degenerate band crossings, respectively. Reproduced with permission of the authors. 




Figure 11: Experimental LDOS for energies close to valence mDPs (a), to the main Dirac point (b) and to conduction mDPs 
(c) in a sample with a 13.4 nm period. The scale bars in all images are 10 nm long, (d): Measured dl/dV (proportional to the 
DOS) with arrows pointing to the dips respective to mDPs. Red (black) corresponds to a 13.4 nm (9.0 nm) Moire wavelength. 
Plots were extracted from Ref. [376] with permission of the authors. 


(iii) a triple degenerate band crossing. Examples of band structures where (i)-(iii) are encountered appear 
in Fig. 10. 

In agreement with these calculations, mDPs were observed in Ref. [376] at energies between 0.2 and 1 eV 
by means of STM measurements. Dips in dl/dV curves like the ones reproduced in Fig. 11 are their main 
signature. The predicted electron hole asymmetry comes from the interaction with the substrate, although 
next-nearest neighbour hopping could also contribute to it. Further experimental evidence of mDPs could 
be extracted from (far) infrared spectra, where dissimilarities with respect to suspended graphene’s constant 
conductivity would show up; or by Hall coefficient measurements, which would vary non monotonically upon 
doping with electrons or holes. 

As for the dispersion relation pinned to mDPs described in (i), further analysis shows that it is anisotropic 
with a new Fermi velocity [376] 

vp(9) = yj {vp cos (50) 2 + [Vsin 56 / (2hGj)] 2 . (101) 

Here, 56 is the angle between Gy/2 and 5k = k — Gj/ 2, and V is the first star amplitude of the adhesion 
potential Fourier expansion. 

Bilayer graphene over hBN was featured in Ref. [374] using a similar approach. mDPs also appear 
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in this system, although the presence of a single cone is rare. As a new feature, gaps between the first 
Moire miniband and the rest of the spectrum are common in this case. The dependence of the spectra on 
the rotation angle is now enhanced with respect to single layer graphene due to higher trigonal warping 
corrections. 

Microscopic derivations of spontaneous deformations .— A step beyond the symmetry-based approaches 
described above was recently taken, [262] and specific values for the spontaneous deformations uq were 
predicted as a function of the relative rotation angle between layers. From this, the specific form of the 
low energy electronic model was also computed. [384] This program requires the knowledge of the adhesion 
potential Vs, whose spatial periodicity is given by that of the substrate, and its amplitude modulation is 
parameterized by adhesion energy differences A cab = £ab — £aa and A cba = £ba — £aa- Here ey denotes 
the adhesion energy per graphene unit cell for local ij stacking. The stacking notation is AA for perfectly 
aligned hexagonal lattices, AB for Carbon-on-Boron and BA for Carbon-on-Nitrogen Bernal stackings. The 
values of the different adhesion energies e,, can be extracted from ab initio calculations, [380, 385, 381, 386] 
which almost invariably predict a stronger adhesion for AR-stacking. Ref. [262] constructs a first star 
adhesion potential corresponding to a given A cba and A cab 


±3 

V s (r f ) = Y Vje^+vo-, 


j=± i 


( 102 ) 


V 3> 0 — v j< 0 


A cab + A cba , .A cab — Ae Bj 4 


18 


6^3 


(103) 


g' belonging to the reciprocal lattice of hBN. The uniform offset Vo = (cab+^ba+^aa )/3 may be disregarded 
in the following. In the continuum limit, the total adhesion energy can be written as 


U s = 




d 2 rV s [r,u(r)\] 


±3 


V s [r,u{rj\ = Y 


Vi e 


-iGj r e ig- ■u(r 


j=±l 


(104) 

(105) 


fig is the area of the superlattice unit cell. The linear distortion regime may be considered, 


V S [r, u{r)} ~ V s [r, 0] + u{r) ■ d a V s [r, u(r)] | S=Q , 


(106) 


This approximation allows us to compute analytically an equilibrium deformation by minimizing the adhesion 
energy U$ plus the elastic energy due to in-plane distortions ^™‘P lane ( see below), and can be shown to be 
legitimate, since the inclusion in Eq. (106) of higher order terms in u(r) does not significantly modify the 
result. [262] 

The elastic energy prevents the layer from conforming completely to Vs- According to continuum 

elasticity theory, the elastic energy due to in-plane stretching, ^7 , ™ _ P lane j f or an isotropic medium is obtained 
from Eq. (10) by replacing the strain tensor u,j by the linear strain tensor upp), 
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Figure 12: TOP: Adhesion energy density Vs [r, u{r)] in real space, relative to the average adhesion uq. The rotation angles 
are 9 = 0 (left), 6 = 1.5° (middle) and 6 = 4° (right). BOTTOM: Theoretical FWHM of a section of the deformation profile 
as a function of the rotation angle or Moire period. 


fl gi analogous to f Is, corresponds to the area of the graphene unit cell. A and p, are the Lame constants. 

The minimization of the total energy of the lattice, U = F s I ( 1 ~ plane with respect to the parameters 

ug yields the equilibrium configuration, 
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( 110 ) 


Fig. (12) plots the resulting adhesion energy Vs[r, u(r)] for different relative rotation angles. It is expected 
to be representative of the Young modulus measured in Ref. [260] . A good qualitative agreement is achieved 
with experiments (see Fig. 6): greater angles give rise to configurations close to an undeformed floating 
phase, whereas a small twist drives the system to sharper deformations without a global expansion or 
contraction (a ‘generalized’ floating phase). The inclusion of e.g. global deformations as additional degrees 
of freedom could be considered, but remain unexplored. Other works have studied also the role of out of plane 
instabilities. [382] These could account for quantitative differences with experimental results, noticeable in 
particular in the normalized width of stacking solitons FWHM/L for long Moire periods. In fact, an abrupt 
jump at L ~ 10 nm was reported therein that is not reproduced by the linear deformation model described 
here. Interestingly, however, some conflicting results exist in this regard, with e.g. Ref. [387] reporting no 
evidence of such discontinuity in electronic observables. 

Microscopic derivations of effective electronic model. - Once an equilibrium deformation solution such as 
Eq. (110) has been established, and bearing also Eqs. (19),(20) in mind, an effective low energy Hamiltonian 
may be derived microscopically for graphene on the crystalline substrate, which has the advantage over 
symmetry-based approaches of predicting specific values of all model parameters. The first step involves 
writing a low energy model including graphene and hBN degrees of freedom, 
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.Hg r denotes the Hamiltonian of monolayer graphene and i?hBN corresponds to a single layer of hBN, that 
is approximated by a perfect insulator. e v (e c ) are the energies of its valence (conduction) band. T(r) is 
the interlayer hopping taking deformations in the linear regime into account, [384] and A I\ joins the Dirac 
point of pristine graphene with the analogous in hBN. hBN orbitals can then be integrated out of Eq. (Ill) 
exactly, which results in a local self energy, 

£(r) ss —T^(r)A[Cg N T'(r) (112) 


that must be added to H gr . Together with the induced pseudopotentials {V s , A el ), the final effective Hamil¬ 
tonian reads [384] 


H cS = H 


gr 


k- jA e \r) 


V s (f) + E(f). 


(113) 


Ref. [384] focuses on unrotated graphene over hBN, explicitly spelling out Eq. (113) after decomposing 
all contributions in terms of cr matrices, and into spatially symmetric and antisymmetric components. Every 
harmonic beyond the first star may be neglected in the model. The solution depends solely on dimensionless 
parameters z u (maximum difference in local relative expansion between different regions, see Ref. [384]), rri± 
(associated to the hBN gap and the interlayer hopping t±), [energy scale of the deformation potential, 
Eq. (19)] and g 2 [relating strains to pseudomagnetic fields, Eq. (20)]. 

The electronic spectrum, density of states (DOS) and local density of states (LDOS) were analysed as a 
function of strains, hopping to hBN, V s and A el . The most notable conclusions are: (i) V s and A el do not 
generate a gap by themselves (i.e., when the interlayer hopping t± —> 0), and in fact tend to cancel each 
other’s contribution when t± ^ 0; (ii) the spectral gap at the primary Dirac point, that can be calculated 
analytically, is linear in the distortions z u to leading order; (iii) new Dirac Points emerge in the valence 
and conduction bands, in agreement with the symmetry-based analysis; (iv) LDOS homogeneity is strongly 
broken by (V Sl A el ); (v) the overall electronic structure at finite energies shows a strong sensitivity to the 
values of {z u ,m±, g±, 52 ) within their realistic range. As an example, the spectral gap can range from 0 to 
~ lOmeV depending on their choice. Fig. 13 illustrates these results. Therefore, a direct measurement of 
electronic structure properties, such as the DOS, could allow to obtain accurate values of (z u , m±, g\, < 72 ), 
to be contrasted with the theoretical predictions. 

Ref. [384] uses the values z u = —0.18 (which corresponds to Ukk /2 = 2.8%, with a graphene lattice 
parameter of a = 2.4 A, lattice mismatch between hBN (ahBN) and graphene of S = OhBN/u — 1 = 0.018, 
A + 2/i = 19.1 eV/A 2 and Acab = —100 meV/unit cell), t± = 0.3 eV and g 2 = 4.8nm _1 . Moreover, t = 3.16 
eV, e c = 3.34 eV and e v = —1.4 eV. As recent works have shown that electronic screening makes graphene’s 
deformation potential cq negligible in practice, [38, 39] g\ = 0 was assumed here. Intriguingly, this choice 
of parameters achieves the best fit to LDOS experimental data[376] when switching off A el but taking into 
account in plane deformations and hoppings to hBN (see Figs, lib and 13d). 


6. Strain in the new families of 2D crystals beyond graphene 

The progresses in the fabrication and characterization of graphene have been successfully applied to the 
isolation of single layers of other 2D materials with properties complementary to graphene [388]. Atomi¬ 
cally thin 2D crystals of boron nitride, transition metal dichalcogenides (TMDs) like M 0 S 2 or WS 2 , black 
phosphorus, etc, have been recently obtained by means of mechanical exfoliation [389, 244, 390], chemical 
vapour deposition (CVD) [391, 392] or laser thinning [393]. Similar to graphene, these other families of 2D 
crystals are extremely strong materials with high elasticity and Young modulus, and they can withstand 
non-hydrostatic (e.g., tensile or shear) stresses up to a significant fraction of their ideal strength without 
inelastic relaxation by plasticity or fracture [394]. Furthermore, the techniques to apply strain in graphene 
can also be applied to other 2D materials, and strain engineering has been proposed as an effective approach 
to continuously tune their electronic and optical properties [395]. In this section, we provide a brief overview 
of the effect of strain in other layered materials different from graphene, with special focus on TMDs and 
black phosphorus. For a comprehensive review on strain engineering in semiconducting 2D crystals, we refer 
the reader to Ref. [66]. 
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Figure 13: Electronic structure given by the model of Eq. (113) for unrotated graphene over hBN. Different combinations of 
physical ingredients have been considered: unstrained graphene on hBN (a), strained graphene on hBN without pseudogauge 
potential (b), strained graphene electronically decoupled from hBN (c), and the complete model of a strained graphene on hBN 
including pseudogauge field (d). Represented are the low-energy band structure, the corresponding total DOS, and the LDOS 
in real space for four energies shown as dashed lines in the band structure and DOS. In the LDOS the colouring corresponds 
to blue for zero, and red for the maximum. We have included the value for the (non-perturbative) gap A at the primary Dirac 
point in each case. 
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Figure 14: (a) Schematic diagram of the fabrication process of wrinkled M 0 S 2 nanolayers. An elastomeric substrate is stretched 
prior depositing M 0 S 2 by mechanical exfoliation. The strain is released afterwards, producing buckling-induced delamination 
of the M 0 S 2 flakes, (b) Optical microscopy image of a wrinkled M 0 S 2 flake fabricated by buckling-induced delamination, (c) 
Atomic force microscopy (AFM) topography image of the region marked by the dashed rectangle in (b). (d) Change in the 
energy of the direct bandgap transition as a function of the strain. The dashed lines show the expected bandgap vs strain 
relationship after accounting for the effect of the finite laser spot size and the funnel effect, (e) Schematic diagram explaining 
the funnel effect due to the nonhomogeneous strain in the wrinkled MoS 2 - (f) Theoretical band structure for nonuniform 
strained M 0 S 2 calculated for a zigzag ribbon. The different panels correspond, from left to right, to the band structure for 
the case of an unstrained ribbon, and for wrinkled ribbons with a 2% and 4 % maximum tensile strain, respectively. In the 
photoluminescence experiments, the wavelength of the emitted light is determined by the direct bandgap transition energy, 
indicated in each panel by a vertical arrow, (g) Colour map of the local density of states as a function of the position along the 
wrinkle profile, calculated theoretically for a wrinkle under non-uniform uniaxial strain. (Figures with permission from Ref. 
[396]) 

6.1. Transition Metal Dichalcogenides 

Among the new families of 2D layered materials, TMDs are recently object of special attention. They are 
semiconductors with a gap in the visible range of the electromagnetic spectrum, present a strong spin-orbit 
coupling, and offer the possibility to control quantum degrees of freedom as the electron spin, and the valley 
and layer pseudospin, among other interesting features[397, 398, 399, 400]. Semiconducting TMDs have 
the form MX 2 , where M = Mo, W is a transition metal and X = S, Se is a clialcogen atom. In their bulk 
configuration, they are composed of two-dimensional X — M — X layers stacked on top of each other, coupled 
by weak van der Waals forces. The metal atoms are ordered in a triangular lattice, each of them bonded to 
six chalcogen atoms located in the top and bottom layers. The electronic band structure of MX 2 changes 
from an indirect band gap for multilayer samples to a direct gap semiconductor for single layers [389]. 

TMDs have a phonon structure highly sensitive to strain. This is revealed by the shift of the Raman 
peaks corresponding to the Ai s out-of-plane mode (where the top and bottom X atoms are moving out of 
plane in opposite directions while M is fixed), and the E 2g in-plane mode (where the M and X atoms are 
moving in-plane in opposite directions) [401]. In particular, applying uniaxial strain lifts the degeneracy of 
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the E^g mode, leading to red shifting and splitting into two distinct peaks for strain of just 1% [396, 402, 403]. 
Strain in the samples can be induced by the application of forces across specific axes of the crystal structure. 
This effect can be used to modify the band structure of TMDs, and in fact a reduction of the band gap 
can be achieved under uniaxial compressive strain across the c-axis of the crystal structure of MX 2 . [404] 
On the other hand, tensile strain is expected to lower the electron effective mass[405] and consequently 
improve electron mobility. A semiconductor to metal transition in single layer M 0 S 2 has been predicted for 
compressive (tensile) biaxial strain of about 15% ( 8 %) [406], which are strengths achievable for these kind 
of 2D crystals [394], First principle calculations have predicted that strain can be induced by depositing 
M 0 S 2 on hexagonal boron nitride, which can lead to a direct-to-indirect gap transition [407]. See Fig. 12 
and Sec. V for a discussion on how superlattices may also induce local variations of the LDOS that is seen 
experimentally by STM. 

Within the framework of the effective mass theory, the effect of strain in the low energy electronic 
spectrum can be easily incorporated by means of group theory. Similarly to the case of graphene, Eq. ( 8 ), 
strain tensor components can be arranged in a scalar [A \) and vectorial {E') irreducible representations 
of D 3h = D 3 ® cr h , the point group of TMDs monolayers of the hexagonal polytype. This vector may be 
interpreted as a pseudo-gauge field that, together with the large spin-orbit coupling provided by M atoms 
and the lack of a centre of inversion in the unit cell (note that D 3 h cannot be decomposed into a direct 
product of the inversion group and the group of rotations of the crystal) can lead to a quantum spin Hall 
effect and time-reversal invariant topological phases in these single-layer materials. For low carrier densities, 
as suggested in Ref. [408], semiconducting TMDs under shear strain will develop spin-polarized Landau 
levels residing in different valleys. For the case of M 0 S 2 , gaps between Landau levels have been estimated 
to scale as hwc/ks ~ 2.1B 0 \T] K, where u c = eB 0 /m* is the cyclotron frequency in terms of the effective 
mass to*, and £?o[T] is the strength of the pseudomagnetic field in Tesla. Considering that, in the case of 
graphene, pseudomagnetic fields of B {) [T] ~ 10 — 10 2 T have been experimentally demonstrated [59], and 
taking into account the intrinsic limitations imposed by the value of the shear modulus and the maximum 
tensile strength of these materials, Landau level gaps of up to « 20 K are within experimental reach. 
Other proposals for different polytypes have appeared in the recent literature [409]. In Ref. [410] it has 
been calculated a low energy k • p Hamiltonian that goes beyond the simple Dirac-like model and captures 
some of the main features of the electronic band dispersion of TMDs under strain, such as the shift of the 
conduction and valence band edges of MX 2 in the presence of homogeneous strain, with the corresponding 
transition from direct to indirect gap, as well as the coupling between spin degrees of freedom and strain. 

Some of the aforementioned theoretical results have been confirmed experimentally, and a change in the 
direct band gap up to ~ 45 meV per 1 % applied strain has been measured [402, 403, 411, 412]. Strain 
engineering has been proposed as a powerful tool to create a broad-band optical funnel in layered 2D 
crystalline semiconductors [413]. Recent experimental observations have shown that, a continuous change 
in strain across an atomically thin sheet of M 0 S 2 , leads to a continuous variation of the optical band gap 
[396]. This could allow not only to capture photons across a wide range of the solar spectrum, but also 
can guide the resulting generated excitons towards the contacts. As sketched in Fig. 14 (a)-(c), buckling 
induced delamination can be used to obtain samples of TMDs with a non-uniform profile of strain, up 
to 2.5% tensile strain. The reduction of the direct band gap with strain, shown in Fig. 14 (d), makes 
that the photo-induced excitons, in order to minimize their energy, move to lower bandgap regions before 
recombining. The nonuniform bandgap profile induced by the local strain of the wrinkles therefore generates 
a trap for photogenerated excitons (funnel effect) with a depth of up to 90 meV, as sketched in Fig. 14 (e). 
The effect of non-uniform strain in the band structure and in the local density of states of M 0 S 2 , shown 
in Fig. 14 (f) and (g), has been theoretically studied with a tight-binding model that properly accounts 
for both, d-orbital contribution from the metal as well as p-orbital contribution from the chalcogen, which 
are relevant for the valence and conduction bands [414, 396, 415]. Recently it has been demonstrated that 
an optoelectronic crystal consisting of a superlattice of artificial atoms can be made from biaxially strained 
single layer M 0 S 2 deposited on a patterned nanocone substrate [416]. 

The large mechanical stretchability and flexibility of single layer TMDs together with the absence of 
inversion symmetry makes them good candidates for high-performance piezoelectric materials. This has 
been recently demonstrated experimentally by Wu et ah, Ref. [417], who have shown that cyclic stretching 
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Figure 15: Single layer BP under uniaxial compression along the out-of-plane direction. Results corresponding to three values 
of h = ho are shown (solid line), along with the electronic band dispersion of the unstrained material (dotted lines). Top and 
side views represent the respective relaxed structure, (d) Change of band gap with height. The original layer thickness used in 
the calculations is 2 /iq . (Reproduced with permission from Ref. [418]) 


and releasing of thin M 0 S 2 flakes with an odd number of atomic layers produces oscillating piezoelectric 
voltage and current outputs, converting mechanical energy into electricity. On the other hand, no output 
is observed for flakes with an even number of layers, which possess inversion symmetry [417]. The strain- 
induced polarization charges in single layer M 0 S 2 can modulate charge carrier transport at a MoS 2 -metal 
barrier and enable enhanced strain sensing. This demonstrates the potential of 2D crystals in powering 
nanodevices and tunable/stretchable electronics/optoelectronics. 

6.2. Black Phosphorus 

Black phosphorus (BP) is another layered material that has been recently synthesized in its single 
layer form, also known as phosphorene [419]. BP is a stable allotrope of phosphorus and an elemental 
semiconductor. At ambient conditions, orthorhombic BP forms puckered layers which are stacked together 
by weak Van der Waals forces (see Fig. 15). Similarly to graphite and TMDs, the layered structure allows 
for mechanical exfoliation to prepare thin layers on a substrate. The band gap is direct, located at the P 
point of the rectangular Brillouin zone, and its size highly depends on the number of layers, suggesting that 
multi-layer BP optical devices could allow photo-detection over a broad spectral range. Furthermore, BP is 
another promising material in nanoelectronics applications since the first field-effect transistors fabricated 
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with few layer samples show drain current modulation on the order of 10 5 (four orders of magnitude larger 
than that in graphene) and high charge-carrier mobilities [420]. 

Before the isolation of single layer BP, the effect of pressure on bulk samples was studied experimentally 
[421] and theoretically [422], The stability of BP nanotubes was studied with density function tight-binding 
methods [423], and a Young modulus of about 300 GPa was estimated. Recent ab initio calculations have 
shown that single layer phosphorene can withstand strain up to 30%, and stress up to 18 GPa and 8 
GPa in the zigzag and armchair directions, respectively [424, 425]. Furthermore, uniaxial stress along the 
direction perpendicular to the layer can be used to change the gap size in the system, eventually driving 
this semiconducting material into a 2D metal [418]. DFT calculations of Ref. [418] considered the effect 
of strain by including in the monolayer unit cell and corresponding atomic positions the constraint z = ±/i 
for all the atoms (see Fig. 15). Therefore an in-plane expansion of the unit cell appears under compressive 
strain (h < ho, where 2/io is the thickness of the free layer). The structure of the BP bonds remains until 
h/ho ~ 0.4, leading to a lattice structure similar to that of a puckered graphene layer. The electronic band 
structure in BP is highly sensitive to strain, suffering a direct to indirect gap transition as we increase the 
compression, and a transition from indirect gap semiconductor to metal is expected for h/h 0 ~ 0.75, as 
shown in Fig. 15. The application of normal compression on bilayer samples of BP can lead to a direct- 
to-indirect gap transition for ~ 3% of strain, and to a fully reversible semiconductor to metal transition at 
~ 13% of applied compression [426]. 

The peculiar puckered structure of BP layers leads to highly anisotropic in-plane optical and electronic 
properties [427, 428]. By applying appropriate uniaxial or biaxial strain in single layer and few-layer BP, 
it has been shown that the anisotropy of the electron effective mass and corresponding mobility direction 
can be rotated by 90° while the anisotropy of holes is not affected by the strain [429] . External strain has 
been shown to drive a Dirac-like dispersion in BP, which differs from the previously reported for graphene, 
silicene, and germanene, in the fact that the anisotropic dispersion of BP allows carriers to behave as 
either massless Dirac or massive carriers, depending on the transport direction along the armchair or zigzag 
axes, respectively [430]. This strain-engineered anisotropic conductance makes BP a promising material for 
mechanical and electronic applications, such as stretchable electrical devices and mechanically controlled 
logical devices. A huge modulation of the gap, up to ~ 1 eV, has been reported in rippled BP samples under 
inhomogeneous strain, suggesting that this material is an excellent candidate to create exciton funnels. In 
this work it has been theoretically predicted that strong confinement of carriers is possible into narrow 
one-dimensional ~ 150 nm wide quantum wires formed along the valleys of the sample ripples, where the 
sample is under compression and the gap is minimum [431]. 

6.3. Other 2D crystals 

Whereas graphene, hBN, TMDs and BP are the most studied 2D crystals due to their demonstrated 
stability, there are other 2D crystals which present interesting features that deserve to be mentioned. Silicene , 
a single layer of Silicon whose atoms are ordered in a low-buckled crystalline structure, is a 2D material 
with non-trivial topological properties. It has been suggested that strain can be used to experimentally 
observe quantum spin Hall effect in this compound [432]. A single layer of Germanium, Germanane, is a 
2D semiconductor with a (irrespective of stacking pattern and thickness) direct band gap and strain can be 
used to drive the membrane into its metallic phase [433]. Two-dimensional layers of Arsenic or arsenene 
have been predicted to be stable in two types of honeycomb structures: buckled and puckered [434], Both 
forms present an indirect gap, that can be tuned into a direct gap by applying strain. Strain induced tuning 
of the band gap has been proven in the transition metal trichalcogenide TiS 3 [435], where it was found that 
single layer and bilayer samples undergo a direct-to-indirect band gap transition when compressive strain is 
induced along the preferred transport axis. 

In summary, we have seen that all the experimental and theoretical techniques developed to study and 
control the effect of strain in graphene are being applied and adapted to other 2D crystals which have 
electronic and optical properties highly sensitive to mechanical deformations. This is broadening a field 
of research which is interesting not only from a fundamental point of view, but also because of potential 
applications in flexible and stretchable electronics and optoelectronics nanodevices. 
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